24 #include "getfem/getfem_generic_assembly_compile_and_exec.h"
25 #include "getfem/getfem_generic_assembly_functions_and_operators.h"
35 void ga_workspace::init() {
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);
42 add_interpolate_transformation
44 add_interpolate_transformation
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;
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);
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;
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);
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));
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));
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));
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))));
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();
120 variables.emplace(name, var_description(
false, &mf, 0,
121 gmm::sub_interval(), &VV, Q));
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)));
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())));
138 bool ga_workspace::is_internal_variable(
const std::string &name)
const {
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)))
145 VAR_SET::const_iterator it = variables.find(name);
146 return it == variables.end() ? false : it->second.is_internal;
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());
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));
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);
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");
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))
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);
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())
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))
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);
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))
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);
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);
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.is_fem_dofs ? it->second.mf : 0;
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);
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);
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.is_fem_dofs ? it->second.mf : 0;
270 const im_data *imd = it->second.imd;
273 return n * mf->get_qdim();
275 return n * imd->tensor_size().total_size();
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);
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.is_fem_dofs ? it->second.mf : 0;
293 const im_data *imd = it->second.imd;
296 bgeot::multi_index mi = mf->get_qdims();
297 if (n > 1 || it->second.qdims.size() > 1) {
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]);
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) {
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]);
317 return it->second.qdims;
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);
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);
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");
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;
359 bool ga_workspace::interpolate_transformation_exists
360 (
const std::string &name)
const {
361 return (md && md->interpolate_transformation_exists(name)) ||
363 parent_workspace->interpolate_transformation_exists(name)) ||
364 (transformations.find(name) != transformations.end());
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);
379 bool ga_workspace::elementary_transformation_exists
380 (
const std::string &name)
const {
381 return (md && md->elementary_transformation_exists(name)) ||
383 parent_workspace->elementary_transformation_exists(name)) ||
384 (elem_transformations.find(name) != elem_transformations.end());
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);
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;
407 bool ga_workspace::secondary_domain_exists
408 (
const std::string &name)
const {
409 return (md && md->secondary_domain_exists(name)) ||
411 parent_workspace->secondary_domain_exists(name)) ||
412 (secondary_domains.find(name) != secondary_domains.end());
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);
429 ga_workspace::register_region(
const mesh &m,
const mesh_region ®ion) {
430 if (&m == &dummy_mesh())
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);
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,
444 bool function_expr, operation_type op_type,
445 const std::string varname_interpolation) {
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))) {
463 if (op_type != ga_workspace::ASSEMBLY)
464 order = add_derivative_order;
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);
476 for (
const ga_workspace::tree_description &td : trees) {
477 if (td.mim == &mim &&
479 td.secondary_domain == tree.secondary_domain &&
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 &&
486 td.operation == op_type &&
487 td.varname_interpolation == varname_interpolation) {
488 ga_tree &ftree = *(td.ptree);
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, ftree.root->children[1]);
494 ga_semantic_analysis(ftree, *
this, m,
495 ref_elt_dim_of_mesh(m,rg),
false, function_expr);
502 ind_tree = trees.size();
504 trees.push_back(tree_description());
505 trees.back().mim = &mim;
507 trees.back().rg = &rg;
508 trees.back().secondary_domain = tree.secondary_domain;
509 trees.back().ptree =
new ga_tree;
510 trees.back().ptree->swap(tree);
511 pga_tree_node root = trees.back().ptree->root;
512 trees.back().name_test1 = root->name_test1;
513 trees.back().name_test2 = root->name_test2;
514 trees.back().interpolate_name_test1 = root->interpolate_name_test1;
515 trees.back().interpolate_name_test2 = root->interpolate_name_test2;
516 trees.back().order = order;
517 trees.back().operation = op_type;
518 trees.back().varname_interpolation = varname_interpolation;
521 if (op_type == ga_workspace::ASSEMBLY && order < add_derivative_order) {
522 std::set<var_trans_pair> expr_variables;
523 ga_extract_variables((remain ? tree : *(trees[ind_tree].ptree)).root,
524 *
this, m, expr_variables,
true);
525 for (
const var_trans_pair &var : expr_variables) {
526 if (!(is_constant(var.varname))) {
527 ga_tree dtree = (remain ? tree : *(trees[ind_tree].ptree));
531 ga_derivative(dtree, *
this, m, var.varname, var.transname, 1+order);
534 ga_semantic_analysis(dtree, *
this, m,
535 ref_elt_dim_of_mesh(m,rg),
false,function_expr);
538 add_tree(dtree, m, mim, rg, expr, add_derivative_order,
539 function_expr, op_type, varname_interpolation);
546 size_type ga_workspace::add_expression(
const std::string &expr,
548 const mesh_region &rg_,
550 const std::string &secondary_dom) {
551 const mesh_region &rg = register_region(mim.linked_mesh(), rg_);
555 std::vector<ga_tree> ltrees(1);
556 ga_read_string(expr, ltrees[0], macro_dictionary());
557 if (secondary_dom.size()) {
558 GMM_ASSERT1(secondary_domain_exists(secondary_dom),
559 "Unknown secondary domain " << secondary_dom);
560 ltrees[0].secondary_domain = secondary_dom;
563 ga_semantic_analysis(ltrees[0], *
this, mim.linked_mesh(),
564 ref_elt_dim_of_mesh(mim.linked_mesh(),rg),
567 GA_TOC(
"First analysis time");
568 if (ltrees[0].root) {
569 if (test1.size() > 1 || test2.size() > 1) {
572 test2.insert(var_trans_pair());
573 ltrees.resize(test1.size()*test2.size(), ltrees[0]);
574 auto ltree = ltrees.begin();
575 for (
const auto &t1 : test1) {
576 for (
const auto &t2 : test2) {
578 if (ntest2 > 0) selected_test2 = t2;
580 ga_semantic_analysis(*ltree, *
this, mim.linked_mesh(),
581 ref_elt_dim_of_mesh(mim.linked_mesh(),rg),
584 if (ltree != ltrees.end()) ++ltree;
587 if (ntest2 == 0) test2.clear();
590 for (ga_tree <ree : ltrees) {
593 max_order = std::max(ltree.root->nb_test_functions(), max_order);
594 add_tree(ltree, mim.linked_mesh(), mim, rg, expr,
595 add_derivative_order,
true);
599 GA_TOC(
"Time for add expression");
603 void ga_workspace::add_function_expression(
const std::string &expr) {
605 ga_read_string(expr, tree, macro_dictionary());
606 ga_semantic_analysis(tree, *
this, dummy_mesh(), 1,
false,
true);
615 void ga_workspace::add_interpolation_expression(
const std::string &expr,
617 const mesh_region &rg_) {
618 const mesh_region &rg = register_region(m, rg_);
620 ga_read_string(expr, tree, macro_dictionary());
621 ga_semantic_analysis(tree, *
this, m, ref_elt_dim_of_mesh(m,rg),
627 ga_workspace::PRE_ASSIGNMENT);
631 void ga_workspace::add_interpolation_expression(
const std::string &expr,
633 const mesh_region &rg_) {
634 const mesh &m = mim.linked_mesh();
635 const mesh_region &rg = register_region(m, rg_);
637 ga_read_string(expr, tree, macro_dictionary());
638 ga_semantic_analysis(tree, *
this, m, ref_elt_dim_of_mesh(m,rg),
641 GMM_ASSERT1(tree.root->nb_test_functions() == 0,
642 "Invalid expression containing test functions");
643 add_tree(tree, m, mim, rg, expr, 0,
false,
644 ga_workspace::PRE_ASSIGNMENT);
648 void ga_workspace::add_assignment_expression
649 (
const std::string &varname,
const std::string &expr,
const mesh_region &rg_,
651 const im_data *imd = associated_im_data(varname);
652 GMM_ASSERT1(imd != 0,
"Only applicable to im_data");
653 const mesh_im &mim = imd->linked_mesh_im();
654 const mesh &m = mim.linked_mesh();
655 const mesh_region &rg = register_region(m, rg_);
657 ga_read_string(expr, tree, macro_dictionary());
658 ga_semantic_analysis(tree, *
this, m, ref_elt_dim_of_mesh(m,rg),
false,
false);
660 GMM_ASSERT1(tree.root->nb_test_functions() == 0,
661 "Invalid expression containing test functions");
662 add_tree(tree, m, mim, rg, expr, order,
false,
663 before ? ga_workspace::PRE_ASSIGNMENT
664 : ga_workspace::POST_ASSIGNMENT,
669 size_type ga_workspace::nb_trees()
const {
return trees.size(); }
671 ga_workspace::tree_description &ga_workspace::tree_info(
size_type i)
674 bool ga_workspace::used_variables(std::vector<std::string> &vl,
675 std::vector<std::string> &vl_test1,
676 std::vector<std::string> &vl_test2,
677 std::vector<std::string> &dl,
680 std::set<var_trans_pair> vll, dll;
681 for (
const std::string &v : vl) vll.insert(var_trans_pair(v,
""));
682 for (
const std::string &d : dl) dll.insert(var_trans_pair(d,
""));
684 for (
const ga_workspace::tree_description &td : trees) {
685 std::set<var_trans_pair> dllaux;
686 bool fv = ga_extract_variables(td.ptree->root, *
this, *(td.m),
689 if (td.order == order)
690 for (
const auto &t : dllaux)
696 if (td.order == order) {
697 if (variable_group_exists(td.name_test1)) {
698 for (
const std::string &t : variable_group(td.name_test1))
699 vll.insert(var_trans_pair(t, td.interpolate_name_test1));
701 vll.insert(var_trans_pair(td.name_test1,
702 td.interpolate_name_test1));
705 for (
const std::string &t : vl_test1)
706 if (td.name_test1 == t)
709 vl_test1.push_back(td.name_test1);
713 if (td.order == order) {
714 if (variable_group_exists(td.name_test1)) {
715 for (
const std::string &t : variable_group(td.name_test1))
716 vll.insert(var_trans_pair(t, td.interpolate_name_test1));
718 vll.insert(var_trans_pair(td.name_test1,
719 td.interpolate_name_test1));
721 if (variable_group_exists(td.name_test2)) {
722 for (
const std::string &t : variable_group(td.name_test2))
723 vll.insert(var_trans_pair(t, td.interpolate_name_test2));
725 vll.insert(var_trans_pair(td.name_test2,
726 td.interpolate_name_test2));
729 for (
size_type j = 0; j < vl_test1.size(); ++j)
730 if ((td.name_test1 == vl_test1[j]) &&
731 (td.name_test2 == vl_test2[j]))
734 vl_test1.push_back(td.name_test1);
735 vl_test2.push_back(td.name_test2);
738 if (fv) islin =
false;
743 for (
const auto &var : vll)
744 if (vl.size() == 0 || var.varname != vl.back())
745 vl.push_back(var.varname);
747 for (
const auto &var : dll)
748 if (dl.size() == 0 || var.varname != dl.back())
749 dl.push_back(var.varname);
754 bool ga_workspace::is_linear(
size_type order) {
755 std::vector<std::string> vl, vl_test1, vl_test2, dl;
756 return used_variables(vl, vl_test1, vl_test2, dl, order);
759 void ga_workspace::define_variable_group(
const std::string &group_name,
760 const std::vector<std::string> &nl) {
761 GMM_ASSERT1(!(variable_exists(group_name)),
"The name of a group of "
762 "variables cannot be the same as a variable name");
764 std::set<const mesh *> ms;
765 bool is_data_ =
false;
766 for (
size_type i = 0; i < nl.size(); ++i) {
768 is_data_ = is_constant(nl[i]);
770 GMM_ASSERT1(is_data_ == is_constant(nl[i]),
771 "It is not possible to mix variables and data in a group");
773 GMM_ASSERT1(variable_exists(nl[i]),
774 "All variables in a group have to exist in the model");
775 const mesh_fem *mf = associated_mf(nl[i]);
776 GMM_ASSERT1(mf,
"Variables in a group should be fem variables");
777 GMM_ASSERT1(ms.find(&(mf->linked_mesh())) == ms.end(),
778 "Two variables in a group cannot share the same mesh");
779 ms.insert(&(mf->linked_mesh()));
781 variable_groups[group_name] = nl;
785 const std::string &ga_workspace::variable_in_group
786 (
const std::string &group_name,
const mesh &m)
const {
787 if (variable_group_exists(group_name)) {
788 for (
const std::string &t : variable_group(group_name))
789 if (&(associated_mf(t)->linked_mesh()) == &m)
791 GMM_ASSERT1(
false,
"No variable in this group for the given mesh");
797 void ga_workspace::assembly(
size_type order,
bool condensation) {
799 const ga_workspace *w =
this;
800 while (w->parent_workspace) w = w->parent_workspace;
801 if (w->md) w->md->nb_dof();
804 ga_instruction_set gis;
805 ga_compile(*
this, gis, order, condensation);
806 GA_TOCTIC(
"Compile time");
808 size_type nb_tot_dof = condensation ? nb_prim_dof + nb_intern_dof
819 if (KQJpr.use_count()) {
823 }
else if (condensation)
824 GMM_ASSERT1(gmm::mat_nrows(*KQJpr) == nb_tot_dof &&
825 gmm::mat_ncols(*KQJpr) == nb_prim_dof,
"Wrong sizes");
829 gmm::resize(col_unreduced_K, nb_tot_dof, nb_tmp_dof);
830 gmm::resize(row_unreduced_K, nb_tmp_dof, nb_prim_dof);
831 gmm::resize(row_col_unreduced_K, nb_tmp_dof, nb_tmp_dof);
838 if (order == 1 || (order == 2 && condensation)) {
839 if (order == 2 && condensation) {
840 GMM_ASSERT1(V->size() == nb_tot_dof,
841 "Wrong size of assembled vector in workspace");
843 gmm::copy(*V, cached_V);
844 gmm::fill(*V, scalar_type(0));
845 }
else if (V.use_count()) {
849 GMM_ASSERT1(V->size() == nb_tot_dof,
850 "Wrong size of assembled vector in workspace");
856 GA_TOCTIC(
"Init time");
858 GA_TOCTIC(
"Exec time");
861 MPI_SUM_VECTOR(assemb_t.as_vector());
862 }
else if (order == 1 || (order == 2 && condensation)) {
864 MPI_SUM_VECTOR(unreduced_V);
870 std::set<std::string> vars_vec_done;
871 std::set<std::pair<std::string, std::string> > vars_mat_done;
872 for (
const auto &term : gis.unreduced_terms) {
873 const std::string &name1 = term.first;
874 const std::string &name2 = term.second;
875 const std::vector<std::string>
876 vg1_(1,name1), vg2_(1,name2),
877 &vg1 = variable_group_exists(name1) ? variable_group(name1) : vg1_,
878 &vg2 = variable_group_exists(name2) ? variable_group(name2) : vg2_;
880 for (
const std::string &vname1 : vg1) {
881 const mesh_fem *mf1 = associated_mf(vname1);
882 if (mf1 && mf1->is_reduced() && vars_vec_done.count(vname1) == 0) {
883 gmm::sub_interval uI1 = temporary_interval_of_variable(vname1),
884 I1 = interval_of_variable(vname1);
886 gmm::sub_vector(unreduced_V, uI1),
887 gmm::sub_vector(*V, I1));
888 vars_vec_done.insert(vname1);
892 for (
const std::string &vname1 : vg1) {
893 for (
const std::string &vname2 : vg2) {
894 const mesh_fem *mf1 = associated_mf(vname1),
895 *mf2 = associated_mf(vname2);
896 if (((mf1 && mf1->is_reduced()) || (mf2 && mf2->is_reduced()))
897 && vars_mat_done.count(std::make_pair(vname1,vname2)) == 0) {
899 uI1 = temporary_interval_of_variable(vname1),
900 uI2 = temporary_interval_of_variable(vname2),
901 I1 = interval_of_variable(vname1),
902 I2 = interval_of_variable(vname2);
903 if (mf1 && mf1->is_reduced()) {
904 if (condensation && vars_vec_done.count(vname1) == 0) {
906 gmm::sub_vector(unreduced_V, uI1),
907 gmm::sub_vector(*V, I1));
908 vars_vec_done.insert(vname1);
910 if (mf2 && mf2->is_reduced()) {
911 model_real_sparse_matrix aux(I1.size(), uI2.size());
912 model_real_row_sparse_matrix M(I1.size(), I2.size());
913 gmm::mult(gmm::transposed(mf1->extension_matrix()),
914 gmm::sub_matrix(row_col_unreduced_K, uI1, uI2),
916 gmm::mult(aux, mf2->extension_matrix(), M);
917 gmm::add(M, gmm::sub_matrix(*K, I1, I2));
918 }
else if (I2.first() < nb_prim_dof) {
919 model_real_sparse_matrix M(I1.size(), I2.size());
920 gmm::mult(gmm::transposed(mf1->extension_matrix()),
921 gmm::sub_matrix(row_unreduced_K, uI1, I2), M);
922 gmm::add(M, gmm::sub_matrix(*K, I1, I2));
925 model_real_row_sparse_matrix M(I1.size(), I2.size());
926 gmm::mult(gmm::sub_matrix(col_unreduced_K, I1, uI2),
927 mf2->extension_matrix(), M);
928 if (I1.first() < nb_prim_dof) {
929 GMM_ASSERT1(I1.last() <= nb_prim_dof,
"Internal error");
930 gmm::add(M, gmm::sub_matrix(*K, I1, I2));
932 gmm::add(M, gmm::sub_matrix(*KQJpr, I1, I2));
935 vars_mat_done.insert(std::make_pair(vname1,vname2));
944 void ga_workspace::set_include_empty_int_points(
bool include) {
945 include_empty_int_pts = include;
948 bool ga_workspace::include_empty_int_points()
const {
949 return include_empty_int_pts;
952 void ga_workspace::add_temporary_interval_for_unreduced_variable
953 (
const std::string &name)
955 if (variable_group_exists(name)) {
956 for (
const std::string &v : variable_group(name))
957 add_temporary_interval_for_unreduced_variable(v);
958 }
else if (tmp_var_intervals.count(name) == 0) {
959 const mesh_fem *mf = associated_mf(name);
960 if (mf && mf->is_reduced()) {
962 tmp_var_intervals[name] = gmm::sub_interval(nb_tmp_dof, nd);
968 void ga_workspace::clear_expressions() { trees.clear(); }
970 void ga_workspace::print(std::ostream &str) {
971 for (
size_type i = 0; i < trees.size(); ++i)
972 if (trees[i].ptree->root) {
973 cout <<
"Expression tree " << i <<
" of order " <<
974 trees[i].ptree->root->nb_test_functions() <<
" :" << endl;
975 ga_print_node(trees[i].ptree->root, str);
980 void ga_workspace::tree_description::copy(
const tree_description& td) {
982 operation = td.operation;
983 varname_interpolation = td.varname_interpolation;
984 name_test1 = td.name_test1;
985 name_test2 = td.name_test2;
986 interpolate_name_test1 = td.interpolate_name_test1;
987 interpolate_name_test2 = td.interpolate_name_test2;
992 if (td.ptree) ptree =
new ga_tree(*(td.ptree));
995 ga_workspace::tree_description &ga_workspace::tree_description::operator =
996 (
const ga_workspace::tree_description& td)
997 {
if (ptree)
delete ptree; ptree = 0;
copy(td);
return *
this; }
998 ga_workspace::tree_description::~tree_description()
999 {
if (ptree)
delete ptree; ptree = 0; }
1002 const inherit var_inherit)
1003 : md(&md_), parent_workspace(0),
1004 with_parent_variables(var_inherit == inherit::ENABLED ||
1005 var_inherit == inherit::ALL),
1006 nb_tmp_dof(0), macro_dict(md_.macro_dictionary())
1009 nb_prim_dof = with_parent_variables ? md->nb_primary_dof() : 0;
1010 nb_intern_dof = with_parent_variables ? md->nb_internal_dof() : 0;
1011 if (var_inherit == inherit::ALL) {
1012 model::varnamelist vlmd;
1013 md->variable_list(vlmd);
1014 for (
const auto &varname : vlmd)
1015 if (md->is_disabled_variable(varname)) {
1016 if (md->is_affine_dependent_variable(varname)) {
1017 std::string orgvarname = md->org_variable(varname);
1018 if (reenabled_var_intervals.count(orgvarname) == 0) {
1019 size_type varsize = gmm::vect_size(md->real_variable(orgvarname));
1020 reenabled_var_intervals[orgvarname]
1021 = gmm::sub_interval (nb_prim_dof, varsize);
1022 nb_prim_dof += varsize;
1024 reenabled_var_intervals[varname]
1025 = reenabled_var_intervals[orgvarname];
1026 }
else if (reenabled_var_intervals.count(varname) == 0) {
1027 size_type varsize = gmm::vect_size(md->real_variable(varname));
1028 reenabled_var_intervals[varname]
1029 = gmm::sub_interval(nb_prim_dof, varsize);
1030 nb_prim_dof += varsize;
1034 first_intern_dof = nb_prim_dof;
1036 ga_workspace::ga_workspace(
const ga_workspace &gaw,
1037 const inherit var_inherit)
1038 : md(0), parent_workspace(&gaw),
1039 with_parent_variables(var_inherit == inherit::ENABLED ||
1040 var_inherit == inherit::ALL),
1041 nb_tmp_dof(0), macro_dict(gaw.macro_dictionary())
1044 nb_prim_dof = with_parent_variables ? gaw.nb_primary_dof() : 0;
1045 nb_intern_dof = with_parent_variables ? gaw.nb_internal_dof() : 0;
1046 first_intern_dof = with_parent_variables ? gaw.first_internal_dof() : 0;
1048 ga_workspace::ga_workspace()
1049 : md(0), parent_workspace(0), with_parent_variables(false),
1050 nb_prim_dof(0), nb_intern_dof(0), first_intern_dof(0), nb_tmp_dof(0)
1052 ga_workspace::~ga_workspace() { clear_expressions(); }
1058 std::string ga_workspace::extract_constant_term(
const mesh &m) {
1059 std::string constant_term;
1060 for (
const ga_workspace::tree_description &td : trees) {
1061 if (td.order == 1) {
1062 ga_tree local_tree = *(td.ptree);
1063 if (local_tree.root)
1064 ga_node_extract_constant_term(local_tree, local_tree.root, *
this, m);
1065 if (local_tree.root)
1066 ga_semantic_analysis(local_tree, *
this, m,
1067 ref_elt_dim_of_mesh(m,-1),
false,
false);
1068 if (local_tree.root && local_tree.root->node_type != GA_NODE_ZERO) {
1069 constant_term +=
"-("+ga_tree_to_string(local_tree)+
")";
1073 return constant_term;
1080 std::string ga_workspace::extract_order0_term() {
1082 for (
const ga_workspace::tree_description &td : trees) {
1083 if (td.order == 0) {
1084 ga_tree &local_tree = *(td.ptree);
1086 term +=
"+("+ga_tree_to_string(local_tree)+
")";
1088 term =
"("+ga_tree_to_string(local_tree)+
")";
1099 std::string ga_workspace::extract_order1_term(
const std::string &varname) {
1101 for (
const ga_workspace::tree_description &td : trees) {
1102 if (td.order == 1 && td.name_test1 == varname) {
1103 ga_tree &local_tree = *(td.ptree);
1105 term +=
"+("+ga_tree_to_string(local_tree)+
")";
1107 term =
"("+ga_tree_to_string(local_tree)+
")";
1117 std::string ga_workspace::extract_Neumann_term(
const std::string &varname) {
1119 for (
const ga_workspace::tree_description &td : trees) {
1120 if (td.order == 1 && td.name_test1 == varname) {
1121 ga_tree &local_tree = *(td.ptree);
1122 if (local_tree.root)
1123 ga_extract_Neumann_term(local_tree, varname, *
this,
1124 local_tree.root, result);