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