GetFEM  5.4.3
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 <[email protected]>
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_ELT_K,
162  GA_NODE_INTERPOLATE_ELT_B,
163  GA_NODE_INTERPOLATE_NORMAL,
164  GA_NODE_INTERPOLATE_DERIVATIVE,
165  GA_NODE_ELEMENTARY,
166  GA_NODE_ELEMENTARY_VAL,
167  GA_NODE_ELEMENTARY_GRAD,
168  GA_NODE_ELEMENTARY_HESS,
169  GA_NODE_ELEMENTARY_DIVERG,
170  GA_NODE_ELEMENTARY_VAL_TEST,
171  GA_NODE_ELEMENTARY_GRAD_TEST,
172  GA_NODE_ELEMENTARY_HESS_TEST,
173  GA_NODE_ELEMENTARY_DIVERG_TEST,
174  GA_NODE_SECONDARY_DOMAIN,
175  GA_NODE_SECONDARY_DOMAIN_VAL,
176  GA_NODE_SECONDARY_DOMAIN_GRAD,
177  GA_NODE_SECONDARY_DOMAIN_HESS,
178  GA_NODE_SECONDARY_DOMAIN_DIVERG,
179  GA_NODE_SECONDARY_DOMAIN_VAL_TEST,
180  GA_NODE_SECONDARY_DOMAIN_GRAD_TEST,
181  GA_NODE_SECONDARY_DOMAIN_HESS_TEST,
182  GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST,
183  GA_NODE_SECONDARY_DOMAIN_X,
184  GA_NODE_SECONDARY_DOMAIN_NORMAL,
185  GA_NODE_XFEM_PLUS,
186  GA_NODE_XFEM_PLUS_VAL,
187  GA_NODE_XFEM_PLUS_GRAD,
188  GA_NODE_XFEM_PLUS_HESS,
189  GA_NODE_XFEM_PLUS_DIVERG,
190  GA_NODE_XFEM_PLUS_VAL_TEST,
191  GA_NODE_XFEM_PLUS_GRAD_TEST,
192  GA_NODE_XFEM_PLUS_HESS_TEST,
193  GA_NODE_XFEM_PLUS_DIVERG_TEST,
194  GA_NODE_XFEM_MINUS,
195  GA_NODE_XFEM_MINUS_VAL,
196  GA_NODE_XFEM_MINUS_GRAD,
197  GA_NODE_XFEM_MINUS_HESS,
198  GA_NODE_XFEM_MINUS_DIVERG,
199  GA_NODE_XFEM_MINUS_VAL_TEST,
200  GA_NODE_XFEM_MINUS_GRAD_TEST,
201  GA_NODE_XFEM_MINUS_HESS_TEST,
202  GA_NODE_XFEM_MINUS_DIVERG_TEST,
203  GA_NODE_ZERO};
204 
205  typedef std::shared_ptr<std::string> pstring;
206  // Print error message indicating the position in the assembly string
207  void ga_throw_error_msg(pstring expr, size_type pos,
208  const std::string &msg);
209 
210 # define ga_throw_error(expr, pos, msg) \
211  { std::stringstream ss; ss << msg; \
212  ga_throw_error_msg(expr, pos, ss.str()); \
213  GMM_ASSERT1(false, "Error in assembly string" ); \
214  }
215 
216  // Structure for the tensor associated with a tree node
217  struct assembly_tensor {
218  bool is_copied;
219  int sparsity_; // 0: plain, 1: vectorized base, 2: vectorised grad, ...
220  size_type qdim_; // Dimension of the vectorization for sparsity tensors
221  base_tensor t;
222  assembly_tensor *tensor_copied;
223 
224  const base_tensor &tensor() const
225  { return (is_copied ? tensor_copied->tensor() : t); }
226  base_tensor &tensor()
227  { return (is_copied ? tensor_copied->tensor() : t); }
228 
229  void set_sparsity(int sp, size_type q)
230  { sparsity_ = sp; qdim_ = q; }
231 
232  size_type qdim() const { return is_copied ? tensor_copied->qdim() : qdim_; }
233 
234  int sparsity() const
235  { return is_copied ? tensor_copied->sparsity() : sparsity_; }
236 
237  inline void set_to_original() { is_copied = false; }
238 
239  inline void set_to_copy(assembly_tensor &t_) {
240  is_copied = true; sparsity_ = t_.sparsity_; qdim_ = t_.qdim_;
241  t = t_.tensor(); tensor_copied = &(t_);
242  }
243 
244  inline void adjust_sizes(const bgeot::multi_index &ssizes)
245  { t.adjust_sizes(ssizes); }
246 
247  inline void adjust_sizes()
248  { if (t.sizes().size() || t.size() != 1) t.init(); }
249 
250  inline void adjust_sizes(size_type i)
251  { if (t.sizes().size() != 1 || t.sizes()[0] != i) t.init(i); }
252 
253  inline void adjust_sizes(size_type i, size_type j) {
254  if (t.sizes().size() != 2 || t.sizes()[0] != i || t.sizes()[1] != j)
255  t.init(i, j);
256  }
257 
258  inline void adjust_sizes(size_type i, size_type j, size_type k) {
259  if (t.sizes().size() != 3 || t.sizes()[0] != i || t.sizes()[1] != j
260  || t.sizes()[2] != k)
261  t.init(i, j, k);
262  }
263  inline void adjust_sizes(size_type i, size_type j,
264  size_type k, size_type l) {
265  if (t.sizes().size() != 3 || t.sizes()[0] != i || t.sizes()[1] != j
266  || t.sizes()[2] != k || t.sizes()[3] != l)
267  t.init(i, j, k, l);
268  }
269 
270  void init_scalar_tensor(scalar_type v)
271  { set_to_original(); t.adjust_sizes(); t[0] = v; }
272 
273  void init_vector_tensor(size_type d)
274  { set_to_original(); t.adjust_sizes(d); }
275 
276  void init_matrix_tensor(size_type n, size_type m)
277  { set_to_original(); t.adjust_sizes(n, m); }
278 
279  void init_identity_matrix_tensor(size_type n) {
280  init_matrix_tensor(n, n);
281  auto itw = t.begin();
282  for (size_type i = 0; i < n; ++i)
283  for (size_type j = 0; j < n; ++j)
284  *itw++ = (i == j) ? scalar_type(1) : scalar_type(0);
285  }
286 
287  void init_third_order_tensor(size_type n, size_type m, size_type l)
288  { set_to_original(); t.adjust_sizes(n, m, l); }
289 
290  void init_fourth_order_tensor(size_type n, size_type m,
291  size_type l, size_type k)
292  { set_to_original(); t.adjust_sizes(n, m, l, k); }
293 
294  const bgeot::multi_index &sizes() const { return t.sizes(); }
295 
296  assembly_tensor()
297  : is_copied(false), sparsity_(0), qdim_(0), tensor_copied(0) {}
298  };
299 
300  struct ga_tree_node;
301  typedef ga_tree_node *pga_tree_node;
302 
303 
304  struct ga_tree_node {
305  GA_NODE_TYPE node_type;
306  GA_TOKEN_TYPE op_type;
307  assembly_tensor t;
308  size_type test_function_type; // -1 = undetermined
309  // 0 = no test function,
310  // 1 = first order, 2 = second order,
311  // 3 = both with always first order in first
312  std::string name_test1, name_test2; // variable names corresponding to test
313  // functions when test_function_type > 0.
314  std::string interpolate_name_test1, interpolate_name_test2; // name
315  // of interpolation transformation if any
316  size_type qdim1, qdim2; // Qdims when test_function_type > 0.
317  size_type nbc1, nbc2, nbc3; // For X (nbc1=coordinate number),
318  // macros (nbc1=param number, nbc2,nbc3 type))
319  // and C_MATRIX (nbc1=order).
320  size_type pos; // Position of the first character in string
321  pstring expr; // Original string, for error messages.
322  std::string name; // variable/constant/function/operator name
323  std::string interpolate_name; // For Interpolate : name of transformation
324  std::string interpolate_name_der; // For Interpolate derivative:
325  // name of transformation
326  std::string elementary_name; // For Elementary_transformation :
327  // name of transformation
328  std::string elementary_target;// For Elementary_transformation :
329  // target variable (for its mesh_fem)
330  size_type der1, der2; // For functions and nonlinear operators,
331  // optional derivative or second derivative.
332  bool symmetric_op;
333  pga_tree_node parent; // Parent node
334  std::vector<pga_tree_node> children; // Children nodes
335  scalar_type hash_value; // Hash value to identify nodes.
336  bool marked; // For specific use of some algorithms
337 
338  inline const base_tensor &tensor() const { return t.tensor(); }
339  inline base_tensor &tensor() { return t.tensor(); }
340  int sparsity() const { return t.sparsity(); }
341 
342  inline size_type nb_test_functions() const {
343  if (test_function_type == size_type(-1)) return 0;
344  return test_function_type - (test_function_type >= 2 ? 1 : 0);
345  }
346 
347  inline size_type tensor_order() const
348  { return t.sizes().size() - nb_test_functions(); }
349 
350  inline size_type tensor_test_size() const {
351  size_type st = nb_test_functions();
352  return (st >= 1 ? t.sizes()[0] : 1) * (st == 2 ? t.sizes()[1] : 1);
353  }
354 
355  inline size_type tensor_proper_size() const
356  { return t.tensor().size() / tensor_test_size(); }
357 
358  inline size_type tensor_proper_size(size_type i) const
359  { return t.sizes()[nb_test_functions()+i]; }
360 
361 
362  void mult_test(const pga_tree_node n0, const pga_tree_node n1);
363 
364  bool tensor_is_zero() {
365  if (node_type == GA_NODE_ZERO) return true;
366  if (node_type != GA_NODE_CONSTANT) return false;
367  for (size_type i = 0; i < tensor().size(); ++i)
368  if (tensor()[i] != scalar_type(0)) return false;
369  return true;
370  }
371 
372  inline bool is_constant() {
373  return (node_type == GA_NODE_CONSTANT ||
374  (node_type == GA_NODE_ZERO && test_function_type == 0));
375  }
376 
377  inline void init_scalar_tensor(scalar_type v)
378  { t.init_scalar_tensor(v); test_function_type = 0; }
379 
380  inline void init_vector_tensor(size_type d)
381  { t.init_vector_tensor(d); test_function_type = 0; }
382 
383  inline void init_matrix_tensor(size_type n, size_type m)
384  { t.init_matrix_tensor(n, m); test_function_type = 0; }
385 
386  inline void init_identity_matrix_tensor(size_type n)
387  { t.init_identity_matrix_tensor(n); test_function_type = 0; }
388 
389  inline void init_third_order_tensor(size_type n, size_type m, size_type l)
390  { t.init_third_order_tensor(n, m, l); test_function_type = 0; }
391 
392  inline void init_fourth_order_tensor(size_type n, size_type m,
393  size_type l, size_type k)
394  { t.init_fourth_order_tensor(n, m, l, k); test_function_type = 0; }
395 
396  inline void adopt_child(pga_tree_node new_child)
397  { children.push_back(new_child); children.back()->parent = this; }
398 
399  inline void accept_child(size_type i) {
400  GMM_ASSERT1(i < children.size(), "Internal error");
401  children[i]->parent = this;
402  }
403 
404  inline void replace_child(pga_tree_node oldchild,
405  pga_tree_node newchild) {
406  bool found = false;
407  for (pga_tree_node &child : children)
408  if (child == oldchild) { child = newchild; found = true; }
409  GMM_ASSERT1(found, "Internal error");
410  }
411 
412  ga_tree_node()
413  : node_type(GA_NODE_VOID), test_function_type(-1), qdim1(0), qdim2(0),
414  nbc1(0), nbc2(0), nbc3(0), pos(0), expr(0), der1(0), der2(0),
415  symmetric_op(false), hash_value(0) {}
416  ga_tree_node(GA_NODE_TYPE ty, size_type p, pstring expr_)
417  : node_type(ty), test_function_type(-1),
418  qdim1(0), qdim2(0), nbc1(0), nbc2(0), nbc3(0),
419  pos(p), expr(expr_), der1(0), der2(0), symmetric_op(false),
420  hash_value(0) {}
421  ga_tree_node(scalar_type v, size_type p, pstring expr_)
422  : node_type(GA_NODE_CONSTANT), test_function_type(-1),
423  qdim1(0), qdim2(0), nbc1(0), nbc2(0), nbc3(0),
424  pos(p), expr(expr_), der1(0), der2(0), symmetric_op(false),
425  hash_value(0)
426  { init_scalar_tensor(v); }
427  ga_tree_node(const char *n, size_type l, size_type p, pstring expr_)
428  : node_type(GA_NODE_NAME), test_function_type(-1),
429  qdim1(0), qdim2(0), nbc1(0), nbc2(0), nbc3(0),
430  pos(p), expr(expr_), name(n, l), der1(0), der2(0), symmetric_op(false),
431  hash_value(0) {}
432  ga_tree_node(GA_TOKEN_TYPE op, size_type p, pstring expr_)
433  : node_type(GA_NODE_OP), op_type(op), test_function_type(-1),
434  qdim1(0), qdim2(0), nbc1(0), nbc2(0), nbc3(0),
435  pos(p), expr(expr_), der1(0), der2(0), symmetric_op(false),
436  hash_value(0) {}
437  };
438 
439  struct ga_tree {
440  pga_tree_node root, current_node;
441  std::string secondary_domain;
442 
443  void add_scalar(scalar_type val, size_type pos, pstring expr);
444  void add_allindices(size_type pos, pstring expr);
445  void add_name(const char *name, size_type length, size_type pos,
446  pstring expr);
447  void add_sub_tree(ga_tree &sub_tree);
448  void add_params(size_type pos, pstring expr);
449  void add_matrix(size_type pos, pstring expr);
450  void add_op(GA_TOKEN_TYPE op_type, size_type pos, pstring expr);
451  void clear_node_rec(pga_tree_node pnode);
452  void clear_node(pga_tree_node pnode);
453  void clear() { clear_node_rec(root); root = current_node = nullptr; }
454  void clear_children(pga_tree_node pnode);
455  void replace_node_by_child(pga_tree_node pnode, size_type i);
456  void copy_node(pga_tree_node pnode, pga_tree_node &pnode_new);
457  void duplicate_with_operation(pga_tree_node pnode, GA_TOKEN_TYPE op_type);
458  void duplicate_with_addition(pga_tree_node pnode)
459  { duplicate_with_operation(pnode, GA_PLUS); }
460  void duplicate_with_substraction(pga_tree_node pnode)
461  { duplicate_with_operation(pnode, GA_MINUS); }
462  void insert_node(pga_tree_node pnode, GA_NODE_TYPE node_type);
463  void add_child(pga_tree_node pnode, GA_NODE_TYPE node_type = GA_NODE_VOID);
464  void swap(ga_tree &tree) {
465  std::swap(root, tree.root);
466  std::swap(current_node, tree.current_node);
467  std::swap(secondary_domain, tree.secondary_domain);
468  }
469 
470  ga_tree() : root(nullptr), current_node(nullptr), secondary_domain() {}
471 
472  ga_tree(const ga_tree &tree)
473  : root(nullptr), current_node(nullptr),
474  secondary_domain(tree.secondary_domain)
475  {
476  if (tree.root)
477  copy_node(tree.root, root);
478  }
479 
480  ga_tree &operator =(const ga_tree &tree) {
481  clear();
482  secondary_domain = tree.secondary_domain;
483  if (tree.root)
484  copy_node(tree.root, root);
485  return *this;
486  }
487 
488  ~ga_tree() { clear(); }
489  };
490 
491  // Test equality or equivalence of two sub trees.
492  // version = 0 : strict equality
493  // 1 : give the same result
494  // 2 : give the same result with transposition of test functions
495  bool sub_tree_are_equal
496  (const pga_tree_node pnode1, const pga_tree_node pnode2,
497  const ga_workspace &workspace, int version);
498 
499  // Transform the expression of a node and its sub-nodes in the equivalent
500  // assembly string sent to ostream str
501  void ga_print_node(const pga_tree_node pnode, std::ostream &str);
502  // The same for the whole tree, the result is a std::string
503  std::string ga_tree_to_string(const ga_tree &tree);
504 
505  // Syntax analysis of an assembly string. Conversion to a tree.
506  // No semantic analysis is done. The tree can be inconsistent.
507  void ga_read_string(const std::string &expr, ga_tree &tree,
508  const ga_macro_dictionary &macro_dict);
509  void ga_read_string_reg(const std::string &expr, ga_tree &tree,
510  ga_macro_dictionary &macro_dict);
511 
512 
513 } /* end of namespace */
514 
515 
516 #endif /* GETFEM_GENERIC_ASSEMBLY_TREE_H__ */
bgeot_geotrans_inv.h
Inversion of geometric transformations.
bgeot_rtree.h
region-tree for window/point search on a set of rectangles.
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_assembly.h
A language for generic assembly of pde boundary value problems.
getfem_models.h
Model representation in Getfem.
getfem_omp.h
Tools for multithreaded, OpenMP and Boost based parallelization.
gmm_blas.h
Basic linear algebra functions.
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
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_accumulated_distro.h:47