GetFEM  5.5
getfem_generic_assembly_workspace.cc
1 /*===========================================================================
2 
3  Copyright (C) 2013-2026 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, see https://www.gnu.org/licenses/.
18 
19 ===========================================================================*/
20 
23 #include "getfem/getfem_generic_assembly_compile_and_exec.h"
24 #include "getfem/getfem_generic_assembly_functions_and_operators.h"
25 
26 namespace getfem {
27 
28  //=========================================================================
29  // Structure dealing with user defined environment : constant, variables,
30  // functions, operators.
31  //=========================================================================
32 
33 
34  void ga_workspace::init() {
35  // allocate own storage for K an V to be used unless/until external
36  // storage is provided with set_assembled_matrix/vector
37  K = std::make_shared<model_real_sparse_matrix>(2,2);
38  V = std::make_shared<base_vector>(2);
39  KQJpr = std::make_shared<model_real_sparse_matrix>(2,2);
40  // add default transformations
41  add_interpolate_transformation // deprecated version
43  add_interpolate_transformation
44  ("neighbor_element", interpolate_transformation_neighbor_instance());
45 
46  ga_tree tree1;
47  pstring s1 = std::make_shared<std::string>("Hess_u");
48  tree1.add_name(s1->c_str(), 6, 0, s1);
49  tree1.root->name = "u";
50  tree1.root->op_type = GA_NAME;
51  tree1.root->node_type = GA_NODE_MACRO_PARAM;
52  tree1.root->nbc1 = 0;
53  tree1.root->nbc2 = ga_parse_prefix_operator(*s1);
54  tree1.root->nbc3 = ga_parse_prefix_test(*s1);
55  ga_macro gam1("Hess", tree1, 1);
56  macro_dict.add_macro(gam1);
57 
58  ga_tree tree2;
59  pstring s2 = std::make_shared<std::string>("Div_u");
60  tree2.add_name(s2->c_str(), 5, 0, s2);
61  tree2.root->name = "u";
62  tree2.root->op_type = GA_NAME;
63  tree2.root->node_type = GA_NODE_MACRO_PARAM;
64  tree2.root->nbc1 = 0;
65  tree2.root->nbc2 = ga_parse_prefix_operator(*s2);
66  tree2.root->nbc3 = ga_parse_prefix_test(*s2);
67  ga_macro gam2("Div", tree2, 1);
68  macro_dict.add_macro(gam2);
69  }
70 
71  // variables and variable groups
72  void ga_workspace::add_fem_variable
73  (const std::string &name, const mesh_fem &mf,
74  const gmm::sub_interval &I, const model_real_plain_vector &VV) {
75  GMM_ASSERT1(nb_intern_dof == 0 || I.last() < first_intern_dof,
76  "The provided interval overlaps with internal dofs");
77  nb_prim_dof = std::max(nb_prim_dof, I.last());
78  variables.emplace(name, var_description(true, &mf, 0, I, &VV, 1));
79  }
80 
81  void ga_workspace::add_im_variable
82  (const std::string &name, const im_data &imd,
83  const gmm::sub_interval &I, const model_real_plain_vector &VV) {
84  GMM_ASSERT1(nb_intern_dof == 0 || I.last() <= first_intern_dof,
85  "The provided interval overlaps with internal dofs");
86  nb_prim_dof = std::max(nb_prim_dof, I.last());
87  variables.emplace(name, var_description(true, 0, &imd, I, &VV, 1));
88  }
89 
90  void ga_workspace::add_internal_im_variable
91  (const std::string &name, const im_data &imd,
92  const gmm::sub_interval &I, const model_real_plain_vector &VV) {
93  GMM_ASSERT1(I.first() >= nb_prim_dof,
94  "The provided interval overlaps with primary dofs");
95  nb_intern_dof += first_intern_dof - std::min(first_intern_dof, I.first());
96  first_intern_dof = std::min(first_intern_dof, I.first());
97  nb_intern_dof += first_intern_dof + nb_intern_dof
98  - std::min(first_intern_dof + nb_intern_dof, I.last());
99  variables.emplace(name, var_description(true, 0, &imd, I, &VV, 1, true));
100  }
101 
102  void ga_workspace::add_fixed_size_variable
103  (const std::string &name,
104  const gmm::sub_interval &I, const model_real_plain_vector &VV) {
105  GMM_ASSERT1(nb_intern_dof == 0 || I.last() <= first_intern_dof,
106  "The provided interval overlaps with internal dofs");
107  nb_prim_dof = std::max(nb_prim_dof, I.last());
108  variables.emplace(name, var_description(true, 0, 0, I, &VV,
109  dim_type(gmm::vect_size(VV))));
110  }
111 
112  void ga_workspace::add_fem_constant
113  (const std::string &name, const mesh_fem &mf,
114  const model_real_plain_vector &VV) {
115  GMM_ASSERT1(mf.nb_dof(), "The provided mesh_fem of variable" << name
116  << "has zero degrees of freedom");
117  size_type Q = gmm::vect_size(VV)/mf.nb_dof();
118  if (Q == 0) Q = size_type(1);
119  variables.emplace(name, var_description(false, &mf, 0,
120  gmm::sub_interval(), &VV, Q));
121  }
122 
123  void ga_workspace::add_fixed_size_constant
124  (const std::string &name, const model_real_plain_vector &VV) {
125  variables.emplace(name, var_description(false, 0, 0,
126  gmm::sub_interval(), &VV,
127  gmm::vect_size(VV)));
128  }
129 
130  void ga_workspace::add_im_data(const std::string &name, const im_data &imd,
131  const model_real_plain_vector &VV) {
132  variables.emplace(name, var_description
133  (false, 0, &imd, gmm::sub_interval(), &VV,
134  gmm::vect_size(VV)/(imd.nb_filtered_index() * imd.nb_tensor_elem())));
135  }
136 
137  bool ga_workspace::is_internal_variable(const std::string &name) const {
138 
139  if ((md && md->variable_exists(name) && md->is_internal_variable(name)) ||
140  (parent_workspace && parent_workspace->variable_exists(name)
141  && parent_workspace->is_internal_variable(name)))
142  return true;
143  else {
144  VAR_SET::const_iterator it = variables.find(name);
145  return it == variables.end() ? false : it->second.is_internal;
146  }
147  }
148 
149  bool ga_workspace::variable_exists(const std::string &name) const {
150  return (md && md->variable_exists(name)) ||
151  (parent_workspace && parent_workspace->variable_exists(name)) ||
152  (variables.find(name) != variables.end());
153  }
154 
155  bool ga_workspace::variable_group_exists(const std::string &name) const {
156  return (variable_groups.find(name) != variable_groups.end()) ||
157  (md && md->variable_group_exists(name)) ||
158  (parent_workspace && parent_workspace->variable_group_exists(name));
159  }
160 
161  const std::vector<std::string>&
162  ga_workspace::variable_group(const std::string &group_name) const {
163  std::map<std::string, std::vector<std::string> >::const_iterator
164  it = variable_groups.find(group_name);
165  if (it != variable_groups.end())
166  return (variable_groups.find(group_name))->second;
167  if (md && md->variable_group_exists(group_name))
168  return md->variable_group(group_name);
169  if (parent_workspace &&
170  parent_workspace->variable_group_exists(group_name))
171  return parent_workspace->variable_group(group_name);
172  GMM_ASSERT1(false, "Undefined variable group " << group_name);
173  }
174 
175  const std::string&
176  ga_workspace::first_variable_of_group(const std::string &name) const {
177  const std::vector<std::string> &t = variable_group(name);
178  GMM_ASSERT1(t.size(), "Variable group " << name << " is empty");
179  return t[0];
180  }
181 
182  bool ga_workspace::is_constant(const std::string &name) const {
183  VAR_SET::const_iterator it = variables.find(name);
184  if (it != variables.end())
185  return !(it->second.is_variable);
186  else if (variable_group_exists(name))
187  return is_constant(first_variable_of_group(name));
188  else if (reenabled_var_intervals.count(name))
189  return false;
190  else if (md && md->variable_exists(name))
191  return md->is_data(name);
192  else if (parent_workspace && parent_workspace->variable_exists(name))
193  return parent_workspace->is_constant(name);
194  GMM_ASSERT1(false, "Undefined variable " << name);
195  }
196 
197  bool ga_workspace::is_disabled_variable(const std::string &name) const {
198  VAR_SET::const_iterator it = variables.find(name);
199  if (it != variables.end())
200  return false;
201  else if (variable_group_exists(name))
202  return is_disabled_variable(first_variable_of_group(name));
203  else if (reenabled_var_intervals.count(name))
204  return false;
205  else if (md && md->variable_exists(name))
206  return md->is_disabled_variable(name);
207  else if (parent_workspace && parent_workspace->variable_exists(name))
208  return parent_workspace->is_disabled_variable(name);
209  GMM_ASSERT1(false, "Undefined variable " << name);
210  }
211 
212  const scalar_type &
213  ga_workspace::factor_of_variable(const std::string &name) const {
214  static const scalar_type one(1);
215  VAR_SET::const_iterator it = variables.find(name);
216  if (it != variables.end()) return one;
217  if (variable_group_exists(name))
218  return one;
219  if (md && md->variable_exists(name)) return md->factor_of_variable(name);
220  if (parent_workspace && parent_workspace->variable_exists(name))
221  return parent_workspace->factor_of_variable(name);
222  GMM_ASSERT1(false, "Undefined variable " << name);
223  }
224 
225  const gmm::sub_interval &
226  ga_workspace::interval_of_variable(const std::string &name) const {
227  VAR_SET::const_iterator it = variables.find(name);
228  if (it != variables.end()) return it->second.I;
229  const auto it2 = reenabled_var_intervals.find(name);
230  if (it2 != reenabled_var_intervals.end()) return it2->second;
231  if (with_parent_variables && md && md->variable_exists(name))
232  return md->interval_of_variable(name);
233  else if (with_parent_variables &&
234  parent_workspace && parent_workspace->variable_exists(name))
235  return parent_workspace->interval_of_variable(name);
236  GMM_ASSERT1(false, "Undefined variable " << name);
237  }
238 
239  const mesh_fem *
240  ga_workspace::associated_mf(const std::string &name) const {
241  VAR_SET::const_iterator it = variables.find(name);
242  if (it != variables.end())
243  return it->second.mf;
244  if (md && md->variable_exists(name))
245  return md->pmesh_fem_of_variable(name);
246  if (parent_workspace && parent_workspace->variable_exists(name))
247  return parent_workspace->associated_mf(name);
248  if (variable_group_exists(name))
249  return associated_mf(first_variable_of_group(name));
250  GMM_ASSERT1(false, "Undefined variable or group " << name);
251  }
252 
253  const im_data *
254  ga_workspace::associated_im_data(const std::string &name) const {
255  VAR_SET::const_iterator it = variables.find(name);
256  if (it != variables.end()) return it->second.imd;
257  if (md && md->variable_exists(name))
258  return md->pim_data_of_variable(name);
259  if (parent_workspace && parent_workspace->variable_exists(name))
260  return parent_workspace->associated_im_data(name);
261  if (variable_group_exists(name)) return 0;
262  GMM_ASSERT1(false, "Undefined variable " << name);
263  }
264 
265  size_type ga_workspace::qdim(const std::string &name) const {
266  VAR_SET::const_iterator it = variables.find(name);
267  if (it != variables.end()) {
268  const mesh_fem *mf = it->second.mf;
269  const im_data *imd = it->second.imd;
270  size_type n = it->second.qdim();
271  if (mf) {
272  return n * mf->get_qdim();
273  } else if (imd) {
274  return n * imd->tensor_size().total_size();
275  }
276  return n;
277  }
278  if (md && md->variable_exists(name))
279  return md->qdim_of_variable(name);
280  if (parent_workspace && parent_workspace->variable_exists(name))
281  return parent_workspace->qdim(name);
282  if (variable_group_exists(name))
283  return qdim(first_variable_of_group(name));
284  GMM_ASSERT1(false, "Undefined variable or group " << name);
285  }
286 
287  bgeot::multi_index
288  ga_workspace::qdims(const std::string &name) const {
289  VAR_SET::const_iterator it = variables.find(name);
290  if (it != variables.end()) {
291  const mesh_fem *mf = it->second.mf;
292  const im_data *imd = it->second.imd;
293  size_type n = it->second.qdim();
294  if (mf) {
295  bgeot::multi_index mi = mf->get_qdims();
296  if (n > 1 || it->second.qdims.size() > 1) {
297  size_type i = 0;
298  if (mi.back() == 1) { mi.back() *= it->second.qdims[0]; ++i; }
299  for (; i < it->second.qdims.size(); ++i)
300  mi.push_back(it->second.qdims[i]);
301  }
302  return mi;
303  } else if (imd) {
304  bgeot::multi_index mi = imd->tensor_size();
305  size_type q = n / imd->nb_filtered_index();
306  GMM_ASSERT1(q % imd->nb_tensor_elem() == 0,
307  "Invalid mesh im data vector");
308  if (n > 1 || it->second.qdims.size() > 1) {
309  size_type i = 0;
310  if (mi.back() == 1) { mi.back() *= it->second.qdims[0]; ++i; }
311  for (; i < it->second.qdims.size(); ++i)
312  mi.push_back(it->second.qdims[i]);
313  }
314  return mi;
315  }
316  return it->second.qdims;
317  }
318  if (md && md->variable_exists(name))
319  return md->qdims_of_variable(name);
320  if (parent_workspace && parent_workspace->variable_exists(name))
321  return parent_workspace->qdims(name);
322  if (variable_group_exists(name))
323  return qdims(first_variable_of_group(name));
324  GMM_ASSERT1(false, "Undefined variable or group " << name);
325  }
326 
327  const model_real_plain_vector &
328  ga_workspace::value(const std::string &name) const {
329  VAR_SET::const_iterator it = variables.find(name);
330  if (it != variables.end())
331  return *(it->second.V);
332  if (md && md->variable_exists(name))
333  return md->real_variable(name);
334  if (parent_workspace && parent_workspace->variable_exists(name))
335  return parent_workspace->value(name);
336  if (variable_group_exists(name))
337  return value(first_variable_of_group(name));
338  GMM_ASSERT1(false, "Undefined variable or group " << name);
339  }
340 
341  scalar_type ga_workspace::get_time_step() const {
342  if (md) return md->get_time_step();
343  if (parent_workspace) return parent_workspace->get_time_step();
344  GMM_ASSERT1(false, "No time step defined here");
345  }
346 
347  void ga_workspace::add_interpolate_transformation
348  (const std::string &name, pinterpolate_transformation ptrans) {
349  if (secondary_domain_exists(name))
350  GMM_ASSERT1(false, "A secondary domain with the same "
351  "name already exists");
352  if (transformations.find(name) != transformations.end())
353  GMM_ASSERT1(name != "neighbor_element", "neighbor_element is a "
354  "reserved interpolate transformation name");
355  transformations[name] = ptrans;
356  }
357 
358  bool ga_workspace::interpolate_transformation_exists
359  (const std::string &name) const {
360  return (md && md->interpolate_transformation_exists(name)) ||
361  (parent_workspace &&
362  parent_workspace->interpolate_transformation_exists(name)) ||
363  (transformations.find(name) != transformations.end());
364  }
365 
366  pinterpolate_transformation
367  ga_workspace::interpolate_transformation(const std::string &name) const {
368  auto it = transformations.find(name);
369  if (it != transformations.end()) return it->second;
370  if (md && md->interpolate_transformation_exists(name))
371  return md->interpolate_transformation(name);
372  if (parent_workspace &&
373  parent_workspace->interpolate_transformation_exists(name))
374  return parent_workspace->interpolate_transformation(name);
375  GMM_ASSERT1(false, "Inexistent transformation " << name);
376  }
377 
378  bool ga_workspace::elementary_transformation_exists
379  (const std::string &name) const {
380  return (md && md->elementary_transformation_exists(name)) ||
381  (parent_workspace &&
382  parent_workspace->elementary_transformation_exists(name)) ||
383  (elem_transformations.find(name) != elem_transformations.end());
384  }
385 
386  pelementary_transformation
387  ga_workspace::elementary_transformation(const std::string &name) const {
388  auto it = elem_transformations.find(name);
389  if (it != elem_transformations.end()) return it->second;
390  if (md && md->elementary_transformation_exists(name))
391  return md->elementary_transformation(name);
392  if (parent_workspace &&
393  parent_workspace->elementary_transformation_exists(name))
394  return parent_workspace->elementary_transformation(name);
395  GMM_ASSERT1(false, "Inexistent elementary transformation " << name);
396  }
397 
398  void ga_workspace::add_secondary_domain(const std::string &name,
399  psecondary_domain psecdom) {
400  if (interpolate_transformation_exists(name))
401  GMM_ASSERT1(false, "An interpolate transformation with the same "
402  "name already exists");
403  secondary_domains[name] = psecdom;
404  }
405 
406  bool ga_workspace::secondary_domain_exists
407  (const std::string &name) const {
408  return (md && md->secondary_domain_exists(name)) ||
409  (parent_workspace &&
410  parent_workspace->secondary_domain_exists(name)) ||
411  (secondary_domains.find(name) != secondary_domains.end());
412  }
413 
414  psecondary_domain
415  ga_workspace::secondary_domain(const std::string &name) const {
416  auto it = secondary_domains.find(name);
417  if (it != secondary_domains.end()) return it->second;
418  if (md && md->secondary_domain_exists(name))
419  return md->secondary_domain(name);
420  if (parent_workspace &&
421  parent_workspace->secondary_domain_exists(name))
422  return parent_workspace->secondary_domain(name);
423  GMM_ASSERT1(false, "Inexistent secondary domain " << name);
424  }
425 
426 
427  const mesh_region &
428  ga_workspace::register_region(const mesh &m, const mesh_region &region) {
429  if (&m == &dummy_mesh())
430  return dummy_mesh_region();
431 
432  std::list<mesh_region> &lmr = registered_mesh_regions[&m];
433  for (const mesh_region &rg : lmr)
434  if (rg.compare(m, region, m)) return rg;
435  lmr.push_back(region);
436  return lmr.back();
437  }
438 
439  void ga_workspace::add_tree(ga_tree &tree, const mesh &m,
440  const mesh_im &mim, const mesh_region &rg,
441  size_type add_derivative_order,
442  bool function_expr, operation_type op_type,
443  const std::string varname_interpolation) {
444  if (tree.root) {
445  // cout << "add tree with tests functions of " << tree.root->name_test1
446  // << " and " << tree.root->name_test2 << endl;
447  // ga_print_node(tree.root, cout); cout << endl;
448 
449  // Eliminate the term if it corresponds to disabled variables
450  if ((tree.root->test_function_type >= 1 &&
451  is_disabled_variable(tree.root->name_test1)) ||
452  (tree.root->test_function_type >= 2 &&
453  is_disabled_variable(tree.root->name_test2))) {
454  // cout<<"disabling term "; ga_print_node(tree.root, cout); cout<<endl;
455  return;
456  }
457 
458  bool remain = true;
459  size_type order = 0, ind_tree = 0;
460 
461  if (op_type != ga_workspace::ASSEMBLY)
462  order = add_derivative_order;
463  else {
464  switch(tree.root->test_function_type) {
465  case 0: order = 0; break;
466  case 1: order = 1; break;
467  case 3: order = 2; break;
468  default: GMM_ASSERT1(false, "Inconsistent term "
469  << tree.root->test_function_type);
470  }
471  }
472 
473  bool found = false;
474  for (const ga_workspace::tree_description &td : trees) {
475  if (td.mim == &mim &&
476  td.m == &m &&
477  td.secondary_domain == tree.secondary_domain &&
478  td.order == order &&
479  td.name_test1 == tree.root->name_test1 &&
480  td.interpolate_name_test1 == tree.root->interpolate_name_test1 &&
481  td.name_test2 == tree.root->name_test2 &&
482  td.interpolate_name_test2 == tree.root->interpolate_name_test2 &&
483  td.rg == &rg &&
484  td.operation == op_type &&
485  td.varname_interpolation == varname_interpolation) {
486  ga_tree &ftree = *(td.ptree);
487 
488  ftree.insert_node(ftree.root, GA_NODE_OP);
489  ftree.root->op_type = GA_PLUS;
490  ftree.root->children.resize(2, nullptr);
491  ftree.copy_node(tree.root, ftree.root->children[1]);
492  ftree.root->children[1]->parent = ftree.root;
493  ga_semantic_analysis(ftree, *this, m,
494  ref_elt_dim_of_mesh(m,rg), false, function_expr);
495  found = true;
496  break;
497  }
498  }
499 
500  if (!found) {
501  ind_tree = trees.size();
502  remain = false;
503  trees.push_back(tree_description());
504  trees.back().mim = &mim;
505  trees.back().m = &m;
506  trees.back().rg = &rg;
507  trees.back().secondary_domain = tree.secondary_domain;
508  trees.back().ptree = new ga_tree;
509  trees.back().ptree->swap(tree);
510  pga_tree_node root = trees.back().ptree->root;
511  trees.back().name_test1 = root->name_test1;
512  trees.back().name_test2 = root->name_test2;
513  trees.back().interpolate_name_test1 = root->interpolate_name_test1;
514  trees.back().interpolate_name_test2 = root->interpolate_name_test2;
515  trees.back().order = order;
516  trees.back().operation = op_type;
517  trees.back().varname_interpolation = varname_interpolation;
518  }
519 
520  if (op_type == ga_workspace::ASSEMBLY && order < add_derivative_order) {
521  std::set<var_trans_pair> expr_variables;
522  ga_extract_variables((remain ? tree : *(trees[ind_tree].ptree)).root,
523  *this, m, expr_variables, true);
524  for (const var_trans_pair &var : expr_variables) {
525  if (!(is_constant(var.varname))) {
526  ga_tree dtree = (remain ? tree : *(trees[ind_tree].ptree));
527  // cout << "Derivation with respect to " << var.varname << " : "
528  // << var.transname << " of " << ga_tree_to_string(dtree) << endl;
529  // GA_TIC;
530  ga_derivative(dtree, *this, m, var.varname, var.transname, 1+order);
531  // cout << "Result : " << ga_tree_to_string(dtree) << endl;
532  // GA_TOCTIC("Derivative time");
533  ga_semantic_analysis(dtree, *this, m,
534  ref_elt_dim_of_mesh(m,rg),false,function_expr);
535  // GA_TOCTIC("Analysis after Derivative time");
536  // cout << "after analysis " << ga_tree_to_string(dtree) << endl;
537  add_tree(dtree, m, mim, rg, add_derivative_order,
538  function_expr, op_type, varname_interpolation);
539  }
540  }
541  }
542  }
543  }
544 
545  size_type ga_workspace::add_expression(const std::string &expr,
546  const mesh_im &mim,
547  const mesh_region &rg_,
548  size_type add_derivative_order,
549  const std::string &secondary_dom) {
550  const mesh_region &rg = register_region(mim.linked_mesh(), rg_);
551  // cout << "adding expression " << expr << endl;
552  GA_TIC;
553  size_type max_order = 0;
554  std::vector<ga_tree> ltrees(1);
555  ga_read_string(expr, ltrees[0], macro_dictionary());
556  if (secondary_dom.size()) {
557  GMM_ASSERT1(secondary_domain_exists(secondary_dom),
558  "Unknown secondary domain " << secondary_dom);
559  ltrees[0].secondary_domain = secondary_dom;
560  }
561  // cout << "read : " << ga_tree_to_string(ltrees[0]) << endl;
562  ga_semantic_analysis(ltrees[0], *this, mim.linked_mesh(),
563  ref_elt_dim_of_mesh(mim.linked_mesh(),rg),
564  false, false, 1);
565  // cout << "analysed : " << ga_tree_to_string(ltrees[0]) << endl;
566  GA_TOC("First analysis time");
567  if (ltrees[0].root) {
568  if (test1.size() > 1 || test2.size() > 1) {
569  size_type ntest2 = test2.size();
570  if (ntest2 == 0) // temporarily add an element to
571  test2.insert(var_trans_pair()); // allow entering the inner loop
572  ltrees.resize(test1.size()*test2.size(), ltrees[0]);
573  auto ltree = ltrees.begin();
574  for (const auto &t1 : test1) {
575  for (const auto &t2 : test2) {
576  selected_test1 = t1;
577  if (ntest2 > 0) selected_test2 = t2;
578  // cout << "analysis with " << selected_test1.first << endl;
579  ga_semantic_analysis(*ltree, *this, mim.linked_mesh(),
580  ref_elt_dim_of_mesh(mim.linked_mesh(),rg),
581  false, false, 2);
582  // cout <<"split: "<< ga_tree_to_string(*ltree) << endl;
583  if (ltree != ltrees.end()) ++ltree;
584  }
585  }
586  if (ntest2 == 0) test2.clear(); // remove temporarily added element
587  }
588 
589  for (ga_tree &ltree : ltrees) {
590  if (ltree.root) {
591  // cout << "adding tree " << ga_tree_to_string(ltree) << endl;
592  max_order = std::max(ltree.root->nb_test_functions(), max_order);
593  add_tree(ltree, mim.linked_mesh(), mim, rg,
594  add_derivative_order, true);
595  }
596  }
597  }
598  GA_TOC("Time for add expression");
599  return max_order;
600  }
601 
602  void ga_workspace::add_function_expression(const std::string &expr) {
603  ga_tree tree;
604  ga_read_string(expr, tree, macro_dictionary());
605  ga_semantic_analysis(tree, *this, dummy_mesh(), 1, false, true);
606  if (tree.root) {
607  // GMM_ASSERT1(tree.root->nb_test_functions() == 0,
608  // "Invalid function expression");
609  add_tree(tree, dummy_mesh(), dummy_mesh_im(), dummy_mesh_region(),
610  0, true);
611  }
612  }
613 
614  void ga_workspace::add_interpolation_expression(const std::string &expr,
615  const mesh &m,
616  const mesh_region &rg_) {
617  const mesh_region &rg = register_region(m, rg_);
618  ga_tree tree;
619  ga_read_string(expr, tree, macro_dictionary());
620  ga_semantic_analysis(tree, *this, m, ref_elt_dim_of_mesh(m,rg),
621  false, false);
622  if (tree.root) {
623  // GMM_ASSERT1(tree.root->nb_test_functions() == 0,
624  // "Invalid expression containing test functions");
625  add_tree(tree, m, dummy_mesh_im(), rg, 0, false,
626  ga_workspace::PRE_ASSIGNMENT);
627  }
628  }
629 
630  void ga_workspace::add_interpolation_expression(const std::string &expr,
631  const mesh_im &mim,
632  const mesh_region &rg_) {
633  const mesh &m = mim.linked_mesh();
634  const mesh_region &rg = register_region(m, rg_);
635  ga_tree tree;
636  ga_read_string(expr, tree, macro_dictionary());
637  ga_semantic_analysis(tree, *this, m, ref_elt_dim_of_mesh(m,rg),
638  false, false);
639  if (tree.root) {
640  GMM_ASSERT1(tree.root->nb_test_functions() == 0,
641  "Invalid expression containing test functions");
642  add_tree(tree, m, mim, rg, 0, false,
643  ga_workspace::PRE_ASSIGNMENT);
644  }
645  }
646 
647  void ga_workspace::add_assignment_expression
648  (const std::string &varname, const std::string &expr, const mesh_region &rg_,
649  size_type order, bool before) {
650  const im_data *imd = associated_im_data(varname);
651  GMM_ASSERT1(imd != 0, "Only applicable to im_data");
652  const mesh_im &mim = imd->linked_mesh_im();
653  const mesh &m = mim.linked_mesh();
654  const mesh_region &rg = register_region(m, rg_);
655  ga_tree tree;
656  ga_read_string(expr, tree, macro_dictionary());
657  ga_semantic_analysis(tree, *this, m, ref_elt_dim_of_mesh(m,rg),false,false);
658  if (tree.root) {
659  GMM_ASSERT1(tree.root->nb_test_functions() == 0,
660  "Invalid expression containing test functions");
661  add_tree(tree, m, mim, rg, order, false,
662  before ? ga_workspace::PRE_ASSIGNMENT
663  : ga_workspace::POST_ASSIGNMENT,
664  varname);
665  }
666  }
667 
668  bool ga_workspace::used_variables(std::vector<std::string> &vl,
669  std::vector<std::string> &vl_test1,
670  std::vector<std::string> &vl_test2,
671  std::vector<std::string> &dl,
672  size_type order) {
673  bool islin = true;
674  std::set<var_trans_pair> vll, dll;
675  for (const std::string &v : vl) vll.insert(var_trans_pair(v, ""));
676  for (const std::string &d : dl) dll.insert(var_trans_pair(d, ""));
677 
678  for (const ga_workspace::tree_description &td : trees) {
679  std::set<var_trans_pair> dllaux;
680  bool fv = ga_extract_variables(td.ptree->root, *this, *(td.m),
681  dllaux, false);
682 
683  if (td.order == order)
684  for (const auto &t : dllaux)
685  dll.insert(t);
686 
687  switch (td.order) {
688  case 0: break;
689  case 1:
690  if (td.order == order) {
691  if (variable_group_exists(td.name_test1)) {
692  for (const std::string &t : variable_group(td.name_test1))
693  vll.insert(var_trans_pair(t, td.interpolate_name_test1));
694  } else {
695  vll.insert(var_trans_pair(td.name_test1,
696  td.interpolate_name_test1));
697  }
698  bool found = false;
699  for (const std::string &t : vl_test1)
700  if (td.name_test1 == t)
701  found = true;
702  if (!found)
703  vl_test1.push_back(td.name_test1);
704  }
705  break;
706  case 2:
707  if (td.order == order) {
708  if (variable_group_exists(td.name_test1)) {
709  for (const std::string &t : variable_group(td.name_test1))
710  vll.insert(var_trans_pair(t, td.interpolate_name_test1));
711  } else {
712  vll.insert(var_trans_pair(td.name_test1,
713  td.interpolate_name_test1));
714  }
715  if (variable_group_exists(td.name_test2)) {
716  for (const std::string &t : variable_group(td.name_test2))
717  vll.insert(var_trans_pair(t, td.interpolate_name_test2));
718  } else {
719  vll.insert(var_trans_pair(td.name_test2,
720  td.interpolate_name_test2));
721  }
722  bool found = false;
723  for (size_type j = 0; j < vl_test1.size(); ++j)
724  if ((td.name_test1 == vl_test1[j]) &&
725  (td.name_test2 == vl_test2[j]))
726  found = true;
727  if (!found) {
728  vl_test1.push_back(td.name_test1);
729  vl_test2.push_back(td.name_test2);
730  }
731  }
732  if (fv) islin = false;
733  break;
734  }
735  }
736  vl.clear();
737  for (const auto &var : vll)
738  if (vl.size() == 0 || var.varname != vl.back())
739  vl.push_back(var.varname);
740  dl.clear();
741  for (const auto &var : dll)
742  if (dl.size() == 0 || var.varname != dl.back())
743  dl.push_back(var.varname);
744 
745  return islin;
746  }
747 
748  bool ga_workspace::is_linear(size_type order) {
749  std::vector<std::string> vl, vl_test1, vl_test2, dl;
750  return used_variables(vl, vl_test1, vl_test2, dl, order);
751  }
752 
753  void ga_workspace::define_variable_group(const std::string &group_name,
754  const std::vector<std::string> &nl) {
755  GMM_ASSERT1(!(variable_exists(group_name)), "The name of a group of "
756  "variables cannot be the same as a variable name");
757 
758  std::set<const mesh *> ms;
759  bool is_data_ = false;
760  for (size_type i = 0; i < nl.size(); ++i) {
761  if (i == 0)
762  is_data_ = is_constant(nl[i]);
763  else {
764  GMM_ASSERT1(is_data_ == is_constant(nl[i]),
765  "It is not possible to mix variables and data in a group");
766  }
767  GMM_ASSERT1(variable_exists(nl[i]),
768  "All variables in a group have to exist in the model");
769  const mesh_fem *mf = associated_mf(nl[i]);
770  GMM_ASSERT1(mf, "Variables in a group should be fem variables");
771  GMM_ASSERT1(ms.find(&(mf->linked_mesh())) == ms.end(),
772  "Two variables in a group cannot share the same mesh");
773  ms.insert(&(mf->linked_mesh()));
774  }
775  variable_groups[group_name] = nl;
776  }
777 
778 
779  const std::string &ga_workspace::variable_in_group
780  (const std::string &group_name, const mesh &m) const {
781  if (variable_group_exists(group_name)) {
782  for (const std::string &t : variable_group(group_name))
783  if (&(associated_mf(t)->linked_mesh()) == &m)
784  return t;
785  GMM_ASSERT1(false, "No variable in this group for the given mesh");
786  } else
787  return group_name;
788  }
789 
790 
791  void ga_workspace::assembly(size_type order, bool condensation) {
792 
793  const ga_workspace *w = this;
794  while (w->parent_workspace) w = w->parent_workspace;
795  if (w->md) w->md->nb_dof(); // To eventually call actualize_sizes()
796 
797  GA_TIC;
798  ga_instruction_set gis;
799  ga_compile(*this, gis, order, condensation);
800  GA_TOCTIC("Compile time");
801 
802  size_type nb_tot_dof = condensation ? nb_prim_dof + nb_intern_dof
803  : nb_prim_dof;
804  if (order == 2) {
805  if (K.use_count()) {
806  gmm::clear(*K);
807  gmm::resize(*K, nb_prim_dof, nb_prim_dof);
808  } // else
809  // We trust that the caller has provided a matrix large enough for the
810  // terms to be assembled (can actually be smaller than the full matrix)
811  // GMM_ASSERT1(gmm::mat_nrows(*K) == nb_prim_dof &&
812  // gmm::mat_ncols(*K) == nb_prim_dof, "Wrong sizes");
813  if (KQJpr.use_count()) {
814  gmm::clear(*KQJpr);
815  if (condensation)
816  gmm::resize(*KQJpr, nb_tot_dof, nb_prim_dof);
817  } else if (condensation)
818  GMM_ASSERT1(gmm::mat_nrows(*KQJpr) == nb_tot_dof &&
819  gmm::mat_ncols(*KQJpr) == nb_prim_dof, "Wrong sizes");
820  gmm::clear(col_unreduced_K);
821  gmm::clear(row_unreduced_K);
822  gmm::clear(row_col_unreduced_K);
823  gmm::resize(col_unreduced_K, nb_tot_dof, nb_tmp_dof);
824  gmm::resize(row_unreduced_K, nb_tmp_dof, nb_prim_dof);
825  gmm::resize(row_col_unreduced_K, nb_tmp_dof, nb_tmp_dof);
826  if (condensation) {
827  gmm::clear(unreduced_V);
828  gmm::resize(unreduced_V, nb_tmp_dof);
829  }
830  }
831 
832  if (order == 1 || (order == 2 && condensation)) {
833  if (order == 2 && condensation) {
834  GMM_ASSERT1(V->size() == nb_tot_dof,
835  "Wrong size of assembled vector in workspace");
836  gmm::resize(cached_V, nb_tot_dof);
837  gmm::copy(*V, cached_V); // current residual is used in condensation
838  gmm::fill(*V, scalar_type(0));
839  } else if (V.use_count()) {
840  gmm::clear(*V);
841  gmm::resize(*V, nb_tot_dof);
842  } else
843  GMM_ASSERT1(V->size() == nb_tot_dof,
844  "Wrong size of assembled vector in workspace");
845  gmm::clear(unreduced_V);
846  gmm::resize(unreduced_V, nb_tmp_dof);
847  }
848  gmm::clear(assembled_tensor().as_vector());
849 
850  GA_TOCTIC("Init time");
851  ga_exec(gis, *this); // --> unreduced_V, *V,
852  GA_TOCTIC("Exec time"); // unreduced_K, *K
853 
854  if (order == 0) {
855  MPI_SUM_VECTOR(assemb_t.as_vector());
856  } else if (order == 1 || (order == 2 && condensation)) {
857  MPI_SUM_VECTOR(*V);
858  MPI_SUM_VECTOR(unreduced_V);
859  }
860 
861  // Deal with reduced fems, unreduced_K --> *K, *KQJpr,
862  // unreduced_V --> *V
863  if (order > 0) {
864  std::set<std::string> vars_vec_done;
865  std::set<std::pair<std::string, std::string> > vars_mat_done;
866  for (const auto &term : gis.unreduced_terms) {
867  const std::string &name1 = term.first;
868  const std::string &name2 = term.second;
869  const std::vector<std::string>
870  vg1_(1,name1), vg2_(1,name2),
871  &vg1 = variable_group_exists(name1) ? variable_group(name1) : vg1_,
872  &vg2 = variable_group_exists(name2) ? variable_group(name2) : vg2_;
873  if (order == 1) {
874  for (const std::string &vname1 : vg1) {
875  const mesh_fem *mf1 = associated_mf(vname1);
876  if (mf1 && mf1->is_reduced() && vars_vec_done.count(vname1) == 0) {
877  gmm::sub_interval uI1 = temporary_interval_of_variable(vname1),
878  I1 = interval_of_variable(vname1);
879  gmm::mult_add(gmm::transposed(mf1->extension_matrix()),
880  gmm::sub_vector(unreduced_V, uI1),
881  gmm::sub_vector(*V, I1));
882  vars_vec_done.insert(vname1);
883  }
884  }
885  } else {
886  for (const std::string &vname1 : vg1) {
887  for (const std::string &vname2 : vg2) {
888  const mesh_fem *mf1 = associated_mf(vname1),
889  *mf2 = associated_mf(vname2);
890  if (((mf1 && mf1->is_reduced()) || (mf2 && mf2->is_reduced()))
891  && vars_mat_done.count(std::make_pair(vname1,vname2)) == 0) {
892  gmm::sub_interval
893  uI1 = temporary_interval_of_variable(vname1),
894  uI2 = temporary_interval_of_variable(vname2),
895  I1 = interval_of_variable(vname1),
896  I2 = interval_of_variable(vname2);
897  if (mf1 && mf1->is_reduced()) {
898  if (condensation && vars_vec_done.count(vname1) == 0) {
899  gmm::mult_add(gmm::transposed(mf1->extension_matrix()),
900  gmm::sub_vector(unreduced_V, uI1),
901  gmm::sub_vector(*V, I1));
902  vars_vec_done.insert(vname1);
903  }
904  if (mf2 && mf2->is_reduced()) {
905  model_real_sparse_matrix aux(I1.size(), uI2.size());
906  model_real_row_sparse_matrix M(I1.size(), I2.size());
907  gmm::mult(gmm::transposed(mf1->extension_matrix()),
908  gmm::sub_matrix(row_col_unreduced_K, uI1, uI2),
909  aux);
910  gmm::mult(aux, mf2->extension_matrix(), M);
911  gmm::add(M, gmm::sub_matrix(*K, I1, I2));
912  } else if (I2.first() < nb_prim_dof) { // !is_internal_variable(vname2)
913  model_real_sparse_matrix M(I1.size(), I2.size());
914  gmm::mult(gmm::transposed(mf1->extension_matrix()),
915  gmm::sub_matrix(row_unreduced_K, uI1, I2), M);
916  gmm::add(M, gmm::sub_matrix(*K, I1, I2));
917  }
918  } else {
919  model_real_row_sparse_matrix M(I1.size(), I2.size());
920  gmm::mult(gmm::sub_matrix(col_unreduced_K, I1, uI2),
921  mf2->extension_matrix(), M);
922  if (I1.first() < nb_prim_dof) {
923  GMM_ASSERT1(I1.last() <= nb_prim_dof, "Internal error");
924  gmm::add(M, gmm::sub_matrix(*K, I1, I2)); // -> *K
925  } else { // vname1 is an internal variable
926  gmm::add(M, gmm::sub_matrix(*KQJpr, I1, I2)); // -> *KQJpr
927  }
928  }
929  vars_mat_done.insert(std::make_pair(vname1,vname2));
930  }
931  }
932  }
933  }
934  }
935  }
936  }
937 
938  void ga_workspace::set_include_empty_int_points(bool include) {
939  include_empty_int_pts = include;
940  }
941 
942  bool ga_workspace::include_empty_int_points() const {
943  return include_empty_int_pts;
944  }
945 
946  void ga_workspace::add_temporary_interval_for_unreduced_variable
947  (const std::string &name)
948  {
949  if (variable_group_exists(name)) {
950  for (const std::string &v : variable_group(name))
951  add_temporary_interval_for_unreduced_variable(v);
952  } else if (tmp_var_intervals.count(name) == 0) {
953  const mesh_fem *mf = associated_mf(name);
954  if (mf && mf->is_reduced()) {
955  size_type nd = mf->nb_basic_dof();
956  tmp_var_intervals[name] = gmm::sub_interval(nb_tmp_dof, nd);
957  nb_tmp_dof += nd;
958  }
959  }
960  }
961 
962  void ga_workspace::clear_expressions() { trees.clear(); }
963 
964  void ga_workspace::print(std::ostream &str) {
965  for (size_type i = 0; i < trees.size(); ++i)
966  if (trees[i].ptree->root) {
967  cout << "Expression tree " << i << " of order " <<
968  trees[i].ptree->root->nb_test_functions() << " :" << endl;
969  ga_print_node(trees[i].ptree->root, str);
970  cout << endl;
971  }
972  }
973 
974  void ga_workspace::tree_description::copy(const tree_description& td) {
975  order = td.order;
976  operation = td.operation;
977  varname_interpolation = td.varname_interpolation;
978  name_test1 = td.name_test1;
979  name_test2 = td.name_test2;
980  interpolate_name_test1 = td.interpolate_name_test1;
981  interpolate_name_test2 = td.interpolate_name_test2;
982  mim = td.mim;
983  m = td.m;
984  rg = td.rg;
985  ptree = 0;
986  if (td.ptree) ptree = new ga_tree(*(td.ptree));
987  }
988 
989  ga_workspace::tree_description &ga_workspace::tree_description::operator =
990  (const ga_workspace::tree_description& td)
991  { if (ptree) delete ptree; ptree = 0; copy(td); return *this; }
992  ga_workspace::tree_description::~tree_description()
993  { if (ptree) delete ptree; ptree = 0; }
994 
995  ga_workspace::ga_workspace(const getfem::model &md_,
996  const inherit var_inherit)
997  : md(&md_), parent_workspace(0),
998  with_parent_variables(var_inherit == inherit::ENABLED ||
999  var_inherit == inherit::ALL),
1000  nb_tmp_dof(0), macro_dict(md_.macro_dictionary())
1001  {
1002  init();
1003  nb_prim_dof = with_parent_variables ? md->nb_primary_dof() : 0;
1004  nb_intern_dof = with_parent_variables ? md->nb_internal_dof() : 0;
1005  if (var_inherit == inherit::ALL) { // enable model's disabled variables
1006  model::varnamelist vlmd;
1007  md->variable_list(vlmd);
1008  for (const auto &varname : vlmd)
1009  if (md->is_disabled_variable(varname)) {
1010  if (md->is_affine_dependent_variable(varname)) {
1011  std::string orgvarname = md->org_variable(varname);
1012  if (reenabled_var_intervals.count(orgvarname) == 0) {
1013  size_type varsize = gmm::vect_size(md->real_variable(orgvarname));
1014  reenabled_var_intervals[orgvarname]
1015  = gmm::sub_interval (nb_prim_dof, varsize);
1016  nb_prim_dof += varsize;
1017  }
1018  reenabled_var_intervals[varname]
1019  = reenabled_var_intervals[orgvarname];
1020  } else if (reenabled_var_intervals.count(varname) == 0) {
1021  size_type varsize = gmm::vect_size(md->real_variable(varname));
1022  reenabled_var_intervals[varname]
1023  = gmm::sub_interval(nb_prim_dof, varsize);
1024  nb_prim_dof += varsize;
1025  }
1026  }
1027  }
1028  first_intern_dof = nb_prim_dof; // dofs are contiguous in getfem::model
1029  }
1030  ga_workspace::ga_workspace(const ga_workspace &gaw,
1031  const inherit var_inherit)
1032  : md(0), parent_workspace(&gaw),
1033  with_parent_variables(var_inherit == inherit::ENABLED ||
1034  var_inherit == inherit::ALL),
1035  nb_tmp_dof(0), macro_dict(gaw.macro_dictionary())
1036  {
1037  init();
1038  nb_prim_dof = with_parent_variables ? gaw.nb_primary_dof() : 0;
1039  nb_intern_dof = with_parent_variables ? gaw.nb_internal_dof() : 0;
1040  first_intern_dof = with_parent_variables ? gaw.first_internal_dof() : 0;
1041  }
1042  ga_workspace::ga_workspace()
1043  : md(0), parent_workspace(0), with_parent_variables(false),
1044  nb_prim_dof(0), nb_intern_dof(0), first_intern_dof(0), nb_tmp_dof(0)
1045  { init(); }
1046  ga_workspace::~ga_workspace() { clear_expressions(); }
1047 
1048  //=========================================================================
1049  // Extract the constant term of degree 1 expressions
1050  //=========================================================================
1051 
1052  std::string ga_workspace::extract_constant_term(const mesh &m) {
1053  std::string constant_term;
1054  for (const ga_workspace::tree_description &td : trees) {
1055  if (td.order == 1) {
1056  ga_tree local_tree = *(td.ptree);
1057  if (local_tree.root)
1058  ga_node_extract_constant_term(local_tree, local_tree.root, *this, m);
1059  if (local_tree.root)
1060  ga_semantic_analysis(local_tree, *this, m,
1061  ref_elt_dim_of_mesh(m,-1), false, false);
1062  if (local_tree.root && local_tree.root->node_type != GA_NODE_ZERO) {
1063  constant_term += "-("+ga_tree_to_string(local_tree)+")";
1064  }
1065  }
1066  }
1067  return constant_term;
1068  }
1069 
1070  //=========================================================================
1071  // Extract the order zero term
1072  //=========================================================================
1073 
1074  std::string ga_workspace::extract_order0_term() {
1075  std::string term;
1076  for (const ga_workspace::tree_description &td : trees) {
1077  if (td.order == 0) {
1078  const ga_tree &local_tree = *(td.ptree);
1079  if (term.size())
1080  term += "+("+ga_tree_to_string(local_tree)+")";
1081  else
1082  term = "("+ga_tree_to_string(local_tree)+")";
1083  }
1084  }
1085  return term;
1086  }
1087 
1088 
1089  //=========================================================================
1090  // Extract the order one term corresponding to a certain test function
1091  //=========================================================================
1092 
1093  std::string ga_workspace::extract_order1_term(const std::string &varname) {
1094  std::string term;
1095  for (const ga_workspace::tree_description &td : trees) {
1096  if (td.order == 1 && td.name_test1 == varname) {
1097  const ga_tree &local_tree = *(td.ptree);
1098  if (term.size())
1099  term += "+("+ga_tree_to_string(local_tree)+")";
1100  else
1101  term = "("+ga_tree_to_string(local_tree)+")";
1102  }
1103  }
1104  return term;
1105  }
1106 
1107  //=========================================================================
1108  // Extract Neumann terms
1109  //=========================================================================
1110 
1111  std::string ga_workspace::extract_Neumann_term(const std::string &varname) {
1112  std::string result;
1113  for (const ga_workspace::tree_description &td : trees) {
1114  if (td.order == 1 && td.name_test1 == varname) {
1115  const ga_tree &local_tree = *(td.ptree);
1116  if (local_tree.root)
1117  ga_extract_Neumann_term(local_tree, varname, *this,
1118  local_tree.root, result);
1119  }
1120  }
1121  return result;
1122  }
1123 
1124 } /* end of namespace */
`‘Model’' variables store the variables, the data and the description of a model.
Semantic analysis of assembly trees and semantic manipulations.
Compilation and execution operations.
void mult_add(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1790
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:976
void fill(L &l, typename gmm::linalg_traits< L >::value_type x)
*‍/
Definition: gmm_blas.h:103
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:58
void resize(V &v, size_type n)
*‍/
Definition: gmm_blas.h:209
void mult(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1663
void add(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:1275
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:48
GEneric Tool for Finite Element Methods.
pinterpolate_transformation interpolate_transformation_neighbor_instance()
Create a new instance of a transformation corresponding to the interpolation on the neighbor element.
const mesh_im & dummy_mesh_im()
Dummy mesh_im for default parameter of functions.
const mesh_region & dummy_mesh_region()
Dummy mesh_region for default parameter of functions.