GetFEM++  5.3
getfem_generic_assembly_tree.h
Go to the documentation of this file.
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  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 
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"
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' operator
88  GA_SKEW, // 'Skew' operator
89  GA_TRACE, // 'Trace' operator
90  GA_DEVIATOR, // 'Deviator' operator
91  GA_INTERPOLATE, // 'Interpolate' operation
92  GA_INTERPOLATE_FILTER, // 'Interpolate_filter' operation
93  GA_ELEMENTARY, // 'Elementary' operation (operation at the element level)
94  GA_XFEM_PLUS, // …valuation on the + side of a level-set for fem_level_set
95  GA_XFEM_MINUS, // …valuation on the - side of a level-set for fem_level_set
96  GA_PRINT, // 'Print' Print the tensor
97  GA_DOT, // '.'
98  GA_DOTMULT, // '.*' componentwise multiplication
99  GA_DOTDIV, // './' componentwise division
100  GA_TMULT, // '@' tensor product
101  GA_COMMA, // ','
102  GA_DCOMMA, // ',,'
103  GA_SEMICOLON, // ';'
104  GA_DSEMICOLON, // ';;'
105  GA_LPAR, // '('
106  GA_RPAR, // ')'
107  GA_LBRACKET, // '['
108  GA_RBRACKET, // ']'
109  GA_NB_TOKEN_TYPE
110  };
111 
112  // Detects Grad_, Hess_ or Div_
113  size_type ga_parse_prefix_operator(std::string &name);
114  // Detects Test_ and Test2_
115  size_type ga_parse_prefix_test(std::string &name);
116 
117  // Types of nodes for the syntax tree
118  enum GA_NODE_TYPE {
119  GA_NODE_VOID = 0,
120  GA_NODE_OP,
121  GA_NODE_PREDEF_FUNC,
122  GA_NODE_SPEC_FUNC,
123  GA_NODE_OPERATOR,
124  GA_NODE_CONSTANT,
125  GA_NODE_NAME,
126  GA_NODE_MACRO_PARAM,
127  GA_NODE_PARAMS,
128  GA_NODE_RESHAPE,
129  GA_NODE_SWAP_IND,
130  GA_NODE_IND_MOVE_LAST,
131  GA_NODE_CONTRACT,
132  GA_NODE_ALLINDICES,
133  GA_NODE_C_MATRIX,
134  GA_NODE_X,
135  GA_NODE_ELT_SIZE,
136  GA_NODE_ELT_K,
137  GA_NODE_ELT_B,
138  GA_NODE_NORMAL,
139  GA_NODE_VAL,
140  GA_NODE_GRAD,
141  GA_NODE_HESS,
142  GA_NODE_DIVERG,
143  GA_NODE_VAL_TEST,
144  GA_NODE_GRAD_TEST,
145  GA_NODE_HESS_TEST,
146  GA_NODE_DIVERG_TEST,
147  GA_NODE_INTERPOLATE,
148  GA_NODE_INTERPOLATE_FILTER,
149  GA_NODE_INTERPOLATE_VAL,
150  GA_NODE_INTERPOLATE_GRAD,
151  GA_NODE_INTERPOLATE_HESS,
152  GA_NODE_INTERPOLATE_DIVERG,
153  GA_NODE_INTERPOLATE_VAL_TEST,
154  GA_NODE_INTERPOLATE_GRAD_TEST,
155  GA_NODE_INTERPOLATE_HESS_TEST,
156  GA_NODE_INTERPOLATE_DIVERG_TEST,
157  GA_NODE_INTERPOLATE_X,
158  GA_NODE_INTERPOLATE_NORMAL,
159  GA_NODE_INTERPOLATE_DERIVATIVE,
160  GA_NODE_ELEMENTARY,
161  GA_NODE_ELEMENTARY_VAL,
162  GA_NODE_ELEMENTARY_GRAD,
163  GA_NODE_ELEMENTARY_HESS,
164  GA_NODE_ELEMENTARY_DIVERG,
165  GA_NODE_ELEMENTARY_VAL_TEST,
166  GA_NODE_ELEMENTARY_GRAD_TEST,
167  GA_NODE_ELEMENTARY_HESS_TEST,
168  GA_NODE_ELEMENTARY_DIVERG_TEST,
169  GA_NODE_XFEM_PLUS,
170  GA_NODE_XFEM_PLUS_VAL,
171  GA_NODE_XFEM_PLUS_GRAD,
172  GA_NODE_XFEM_PLUS_HESS,
173  GA_NODE_XFEM_PLUS_DIVERG,
174  GA_NODE_XFEM_PLUS_VAL_TEST,
175  GA_NODE_XFEM_PLUS_GRAD_TEST,
176  GA_NODE_XFEM_PLUS_HESS_TEST,
177  GA_NODE_XFEM_PLUS_DIVERG_TEST,
178  GA_NODE_XFEM_MINUS,
179  GA_NODE_XFEM_MINUS_VAL,
180  GA_NODE_XFEM_MINUS_GRAD,
181  GA_NODE_XFEM_MINUS_HESS,
182  GA_NODE_XFEM_MINUS_DIVERG,
183  GA_NODE_XFEM_MINUS_VAL_TEST,
184  GA_NODE_XFEM_MINUS_GRAD_TEST,
185  GA_NODE_XFEM_MINUS_HESS_TEST,
186  GA_NODE_XFEM_MINUS_DIVERG_TEST,
187  GA_NODE_ZERO};
188 
189  typedef std::shared_ptr<std::string> pstring;
190  // Print error message indicating the position in the assembly string
191  void ga_throw_error_msg(pstring expr, size_type pos,
192  const std::string &msg);
193 
194 # define ga_throw_error(expr, pos, msg) \
195  { std::stringstream ss; ss << msg; \
196  ga_throw_error_msg(expr, pos, ss.str()); \
197  GMM_ASSERT1(false, "Error in assembly string" ); \
198  }
199 
200  // Structure for the tensor associated with a tree node
201  struct assembly_tensor {
202  bool is_copied;
203  int sparsity_; // 0: plain, 1: vectorized base, 2: vectorised grad, ...
204  size_type qdim_; // Dimension of the vectorization for sparsity tensors
205  base_tensor t;
206  assembly_tensor *tensor_copied;
207 
208  const base_tensor &org_tensor() const
209  { return is_copied ? tensor_copied->org_tensor() : t; }
210  base_tensor &org_tensor()
211  { return is_copied ? tensor_copied->org_tensor() : t; }
212 
213  const base_tensor &tensor() const
214  { return (is_copied ? tensor_copied->tensor() : t); }
215 
216  base_tensor &tensor()
217  { return (is_copied ? tensor_copied->tensor() : t); }
218 
219  void set_sparsity(int sp, size_type q)
220  { sparsity_ = sp; qdim_ = q; }
221 
222  size_type qdim() { return is_copied ? tensor_copied->qdim() : qdim_; }
223 
224  int sparsity() const
225  { return is_copied ? tensor_copied->sparsity() : sparsity_; }
226 
227  inline void set_to_original() { is_copied = false; }
228  inline void set_to_copy(assembly_tensor &t_) {
229  is_copied = true; sparsity_ = t_.sparsity_; qdim_ = t_.qdim_;
230  t = t_.org_tensor(); tensor_copied = &(t_);
231  }
232 
233  inline void adjust_sizes(const bgeot::multi_index &ssizes)
234  { t.adjust_sizes(ssizes); }
235 
236  inline void adjust_sizes()
237  { if (t.sizes().size() || t.size() != 1) t.init(); }
238 
239  inline void adjust_sizes(size_type i)
240  { if (t.sizes().size() != 1 || t.sizes()[0] != i) t.init(i); }
241 
242  inline void adjust_sizes(size_type i, size_type j) {
243  if (t.sizes().size() != 2 || t.sizes()[0] != i || t.sizes()[1] != j)
244  t.init(i, j);
245  }
246 
247  inline void adjust_sizes(size_type i, size_type j, size_type k) {
248  if (t.sizes().size() != 3 || t.sizes()[0] != i || t.sizes()[1] != j
249  || t.sizes()[2] != k)
250  t.init(i, j, k);
251  }
252  inline void adjust_sizes(size_type i, size_type j,
253  size_type k, size_type l) {
254  if (t.sizes().size() != 3 || t.sizes()[0] != i || t.sizes()[1] != j
255  || t.sizes()[2] != k || t.sizes()[3] != l)
256  t.init(i, j, k, l);
257  }
258 
259  void init_scalar_tensor(scalar_type v)
260  { set_to_original(); t.adjust_sizes(); t[0] = v; }
261 
262  void init_vector_tensor(size_type d)
263  { set_to_original(); t.adjust_sizes(d); }
264 
265  void init_matrix_tensor(size_type n, size_type m)
266  { set_to_original(); t.adjust_sizes(n, m); }
267 
268  void init_third_order_tensor(size_type n, size_type m, size_type l)
269  { set_to_original(); t.adjust_sizes(n, m, l); }
270 
271  void init_fourth_order_tensor(size_type n, size_type m,
272  size_type l, size_type k)
273  { set_to_original(); t.adjust_sizes(n, m, l, k); }
274 
275  const bgeot::multi_index &sizes() const { return t.sizes(); }
276 
277  assembly_tensor()
278  : is_copied(false), sparsity_(0), qdim_(0), tensor_copied(0) {}
279  };
280 
281  struct ga_tree_node;
282  typedef ga_tree_node *pga_tree_node;
283 
284 
285  struct ga_tree_node {
286  GA_NODE_TYPE node_type;
287  GA_TOKEN_TYPE op_type;
288  assembly_tensor t;
289  size_type test_function_type; // -1 = undetermined
290  // 0 = no test function,
291  // 1 = first order, 2 = second order,
292  // 3 = both with always first order in first
293  std::string name_test1, name_test2; // variable names corresponding to test
294  // functions when test_function_type > 0.
295  std::string interpolate_name_test1, interpolate_name_test2; // name
296  // of interpolation transformation if any
297  size_type qdim1, qdim2; // Qdims when test_function_type > 0.
298  size_type nbc1, nbc2, nbc3; // For X (nbc1=coordinate number),
299  // macros (nbc1=param number, nbc2,nbc3 type))
300  // and C_MATRIX (nbc1=order).
301  size_type pos; // Position of the first character in string
302  pstring expr; // Original string, for error messages.
303  std::string name; // variable/constant/function/operator name
304  std::string interpolate_name; // For Interpolate : name of transformation
305  std::string interpolate_name_der; // For Interpolate derivative:
306  // name of transformation
307  std::string elementary_name; // For Elementary_transformation :
308  // name of transformation
309  size_type der1, der2; // For functions and nonlinear operators,
310  // optional derivative or second derivative.
311  bool symmetric_op;
312  pga_tree_node parent; // Parent node
313  std::vector<pga_tree_node> children; // Children nodes
314  scalar_type hash_value; // Hash value to identify nodes.
315  bool marked; // For specific use of some algorithms
316 
317  inline const base_tensor &tensor() const { return t.tensor(); }
318  inline base_tensor &tensor() { return t.tensor(); }
319  int sparsity() const { return t.sparsity(); }
320 
321  inline size_type nb_test_functions() const {
322  if (test_function_type == size_type(-1)) return 0;
323  return test_function_type - (test_function_type >= 2 ? 1 : 0);
324  }
325 
326  inline size_type tensor_order() const
327  { return t.sizes().size() - nb_test_functions(); }
328 
329  inline size_type tensor_test_size() const {
330  size_type st = nb_test_functions();
331  return (st >= 1 ? t.sizes()[0] : 1) * (st == 2 ? t.sizes()[1] : 1);
332  }
333 
334  inline size_type tensor_proper_size() const
335  { return t.org_tensor().size() / tensor_test_size(); }
336 
337  inline size_type tensor_proper_size(size_type i) const
338  { return t.sizes()[nb_test_functions()+i]; }
339 
340 
341  void mult_test(const pga_tree_node n0, const pga_tree_node n1);
342 
343  bool tensor_is_zero() {
344  if (node_type == GA_NODE_ZERO) return true;
345  if (node_type != GA_NODE_CONSTANT) return false;
346  for (size_type i = 0; i < tensor().size(); ++i)
347  if (tensor()[i] != scalar_type(0)) return false;
348  return true;
349  }
350 
351  inline void init_scalar_tensor(scalar_type v)
352  { t.init_scalar_tensor(v); test_function_type = 0; }
353 
354  inline void init_vector_tensor(size_type d)
355  { t.init_vector_tensor(d); test_function_type = 0; }
356 
357  inline void init_matrix_tensor(size_type n, size_type m)
358  { t.init_matrix_tensor(n, m); test_function_type = 0; }
359 
360  inline void init_third_order_tensor(size_type n, size_type m, size_type l)
361  { t.init_third_order_tensor(n, m, l); test_function_type = 0; }
362 
363  inline void init_fourth_order_tensor(size_type n, size_type m,
364  size_type l, size_type k)
365  { t.init_fourth_order_tensor(n, m, l, k); test_function_type = 0; }
366 
367  inline void adopt_child(pga_tree_node new_child)
368  { children.push_back(new_child); children.back()->parent = this; }
369 
370  inline void replace_child(pga_tree_node oldchild,
371  pga_tree_node newchild) {
372  bool found = false;
373  for (pga_tree_node &child : children)
374  if (child == oldchild) { child = newchild; found = true; }
375  GMM_ASSERT1(found, "Internal error");
376  }
377 
378  ga_tree_node()
379  : node_type(GA_NODE_VOID), test_function_type(-1), qdim1(0), qdim2(0),
380  nbc1(0), nbc2(0), nbc3(0), pos(0), expr(0), der1(0), der2(0),
381  symmetric_op(false), hash_value(0) {}
382  ga_tree_node(GA_NODE_TYPE ty, size_type p, pstring expr_)
383  : node_type(ty), test_function_type(-1),
384  qdim1(0), qdim2(0), nbc1(0), nbc2(0), nbc3(0),
385  pos(p), expr(expr_), der1(0), der2(0), symmetric_op(false),
386  hash_value(0) {}
387  ga_tree_node(scalar_type v, size_type p, pstring expr_)
388  : node_type(GA_NODE_CONSTANT), test_function_type(-1),
389  qdim1(0), qdim2(0), nbc1(0), nbc2(0), nbc3(0),
390  pos(p), expr(expr_), der1(0), der2(0), symmetric_op(false),
391  hash_value(0)
392  { init_scalar_tensor(v); }
393  ga_tree_node(const char *n, size_type l, size_type p, pstring expr_)
394  : node_type(GA_NODE_NAME), test_function_type(-1),
395  qdim1(0), qdim2(0), nbc1(0), nbc2(0), nbc3(0),
396  pos(p), expr(expr_), name(n, l), der1(0), der2(0), symmetric_op(false),
397  hash_value(0) {}
398  ga_tree_node(GA_TOKEN_TYPE op, size_type p, pstring expr_)
399  : node_type(GA_NODE_OP), op_type(op), test_function_type(-1),
400  qdim1(0), qdim2(0), nbc1(0), nbc2(0), nbc3(0),
401  pos(p), expr(expr_), der1(0), der2(0), symmetric_op(false),
402  hash_value(0) {}
403  };
404 
405  struct ga_tree {
406  pga_tree_node root, current_node;
407 
408  void add_scalar(scalar_type val, size_type pos, pstring expr);
409  void add_allindices(size_type pos, pstring expr);
410  void add_name(const char *name, size_type length, size_type pos,
411  pstring expr);
412  void add_sub_tree(ga_tree &sub_tree);
413  void add_params(size_type pos, pstring expr);
414  void add_matrix(size_type pos, pstring expr);
415  void add_op(GA_TOKEN_TYPE op_type, size_type pos, pstring expr);
416  void clear_node_rec(pga_tree_node pnode);
417  void clear_node(pga_tree_node pnode);
418  void clear() { clear_node_rec(root); root = current_node = nullptr; }
419  void clear_children(pga_tree_node pnode);
420  void replace_node_by_child(pga_tree_node pnode, size_type i);
421  void copy_node(pga_tree_node pnode, pga_tree_node parent,
422  pga_tree_node &child);
423  void duplicate_with_operation(pga_tree_node pnode, GA_TOKEN_TYPE op_type);
424  void duplicate_with_addition(pga_tree_node pnode)
425  { duplicate_with_operation(pnode, GA_PLUS); }
426  void duplicate_with_substraction(pga_tree_node pnode)
427  { duplicate_with_operation(pnode, GA_MINUS); }
428  void insert_node(pga_tree_node pnode, GA_NODE_TYPE node_type);
429  void add_child(pga_tree_node pnode, GA_NODE_TYPE node_type = GA_NODE_VOID);
430  void swap(ga_tree &tree)
431  { std::swap(root, tree.root); std::swap(current_node, tree.current_node); }
432 
433  ga_tree() : root(nullptr), current_node(nullptr) {}
434 
435  ga_tree(const ga_tree &tree) : root(nullptr), current_node(nullptr)
436  { if (tree.root) copy_node(tree.root, nullptr, root); }
437 
438  ga_tree &operator = (const ga_tree &tree)
439  { clear(); if (tree.root) copy_node(tree.root,nullptr,root); return *this; }
440 
441  ~ga_tree() { clear(); }
442  };
443 
444  // Test equality or equivalence of two sub trees.
445  // version = 0 : strict equality
446  // 1 : give the same result
447  // 2 : give the same result with transposition of test functions
448  bool sub_tree_are_equal
449  (const pga_tree_node pnode1, const pga_tree_node pnode2,
450  const ga_workspace &workspace, int version);
451 
452  // Transform the expression of a node and its sub-nodes in the equivalent
453  // assembly string sent to ostream str
454  void ga_print_node(const pga_tree_node pnode,
455  std::ostream &str);
456  // The same for the whole tree, the result is a std::string
457  std::string ga_tree_to_string(const ga_tree &tree);
458 
459  // Syntax analysis of an assembly string. Conversion to a tree.
460  // No semantic analysis is done. The tree can be inconsistent.
461  void ga_read_string(const std::string &expr, ga_tree &tree,
462  const ga_macro_dictionnary &macro_dict);
463  void ga_read_string_reg(const std::string &expr, ga_tree &tree,
464  ga_macro_dictionnary &macro_dict);
465 
466 
467 } /* end of namespace */
468 
469 
470 #endif /* GETFEM_GENERIC_ASSEMBLY_TREE_H__ */
A smart pointer that copies the value it points to on copy operations.
Tools for multithreaded, OpenMP and Boost based parallelization.
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
Inversion of geometric transformations.
Model representation in Getfem.
A simple singleton implementation.
A langage for generic assembly of pde boundary value problems.
GEneric Tool for Finite Element Methods.
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
Basic linear algebra functions.
region-tree for window/point search on a set of rectangles.