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.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);
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.mf;
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.mf;
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->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);
503 ind_tree = trees.size();
505 trees.push_back(tree_description());
506 trees.back().mim = &mim;
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;
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));
532 ga_derivative(dtree, *
this, m, var.varname, var.transname, 1+order);
535 ga_semantic_analysis(dtree, *
this, m,
536 ref_elt_dim_of_mesh(m,rg),
false,function_expr);
539 add_tree(dtree, m, mim, rg, expr, add_derivative_order,
540 function_expr, op_type, varname_interpolation);
547 size_type ga_workspace::add_expression(
const std::string &expr,
549 const mesh_region &rg_,
551 const std::string &secondary_dom) {
552 const mesh_region &rg = register_region(mim.linked_mesh(), rg_);
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;
564 ga_semantic_analysis(ltrees[0], *
this, mim.linked_mesh(),
565 ref_elt_dim_of_mesh(mim.linked_mesh(),rg),
568 GA_TOC(
"First analysis time");
569 if (ltrees[0].root) {
570 if (test1.size() > 1 || test2.size() > 1) {
573 test2.insert(var_trans_pair());
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) {
579 if (ntest2 > 0) selected_test2 = t2;
581 ga_semantic_analysis(*ltree, *
this, mim.linked_mesh(),
582 ref_elt_dim_of_mesh(mim.linked_mesh(),rg),
585 if (ltree != ltrees.end()) ++ltree;
588 if (ntest2 == 0) test2.clear();
591 for (ga_tree <ree : ltrees) {
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);
600 GA_TOC(
"Time for add expression");
604 void ga_workspace::add_function_expression(
const std::string &expr) {
606 ga_read_string(expr, tree, macro_dictionary());
607 ga_semantic_analysis(tree, *
this, dummy_mesh(), 1,
false,
true);
616 void ga_workspace::add_interpolation_expression(
const std::string &expr,
618 const mesh_region &rg_) {
619 const mesh_region &rg = register_region(m, rg_);
621 ga_read_string(expr, tree, macro_dictionary());
622 ga_semantic_analysis(tree, *
this, m, ref_elt_dim_of_mesh(m,rg),
628 ga_workspace::PRE_ASSIGNMENT);
632 void ga_workspace::add_interpolation_expression(
const std::string &expr,
634 const mesh_region &rg_) {
635 const mesh &m = mim.linked_mesh();
636 const mesh_region &rg = register_region(m, rg_);
638 ga_read_string(expr, tree, macro_dictionary());
639 ga_semantic_analysis(tree, *
this, m, ref_elt_dim_of_mesh(m,rg),
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);
649 void ga_workspace::add_assignment_expression
650 (
const std::string &varname,
const std::string &expr,
const mesh_region &rg_,
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_);
658 ga_read_string(expr, tree, macro_dictionary());
659 ga_semantic_analysis(tree, *
this, m, ref_elt_dim_of_mesh(m,rg),
false,
false);
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,
670 size_type ga_workspace::nb_trees()
const {
return trees.size(); }
672 ga_workspace::tree_description &ga_workspace::tree_info(
size_type i)
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,
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,
""));
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),
690 if (td.order == order)
691 for (
const auto &t : dllaux)
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));
702 vll.insert(var_trans_pair(td.name_test1,
703 td.interpolate_name_test1));
706 for (
const std::string &t : vl_test1)
707 if (td.name_test1 == t)
710 vl_test1.push_back(td.name_test1);
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));
719 vll.insert(var_trans_pair(td.name_test1,
720 td.interpolate_name_test1));
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));
726 vll.insert(var_trans_pair(td.name_test2,
727 td.interpolate_name_test2));
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]))
735 vl_test1.push_back(td.name_test1);
736 vl_test2.push_back(td.name_test2);
739 if (fv) islin =
false;
744 for (
const auto &var : vll)
745 if (vl.size() == 0 || var.varname != vl.back())
746 vl.push_back(var.varname);
748 for (
const auto &var : dll)
749 if (dl.size() == 0 || var.varname != dl.back())
750 dl.push_back(var.varname);
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);
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");
765 std::set<const mesh *> ms;
766 bool is_data_ =
false;
767 for (
size_type i = 0; i < nl.size(); ++i) {
769 is_data_ = is_constant(nl[i]);
771 GMM_ASSERT1(is_data_ == is_constant(nl[i]),
772 "It is not possible to mix variables and data in a group");
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()));
782 variable_groups[group_name] = nl;
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)
792 GMM_ASSERT1(
false,
"No variable in this group for the given mesh");
798 void ga_workspace::assembly(
size_type order,
bool condensation) {
800 const ga_workspace *w =
this;
801 while (w->parent_workspace) w = w->parent_workspace;
802 if (w->md) w->md->nb_dof();
805 ga_instruction_set gis;
806 ga_compile(*
this, gis, order, condensation);
807 GA_TOCTIC(
"Compile time");
809 size_type nb_tot_dof = condensation ? nb_prim_dof + nb_intern_dof
820 if (KQJpr.use_count()) {
824 }
else if (condensation)
825 GMM_ASSERT1(gmm::mat_nrows(*KQJpr) == nb_tot_dof &&
826 gmm::mat_ncols(*KQJpr) == nb_prim_dof,
"Wrong sizes");
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);
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");
846 }
else if (V.use_count()) {
850 GMM_ASSERT1(V->size() == nb_tot_dof,
851 "Wrong size of assembled vector in workspace");
857 GA_TOCTIC(
"Init time");
859 GA_TOCTIC(
"Exec time");
862 MPI_SUM_VECTOR(assemb_t.as_vector());
863 }
else if (order == 1 || (order == 2 && condensation)) {
865 MPI_SUM_VECTOR(unreduced_V);
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_;
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);
887 gmm::sub_vector(unreduced_V, uI1),
888 gmm::sub_vector(*V, I1));
889 vars_vec_done.insert(vname1);
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) {
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) {
907 gmm::sub_vector(unreduced_V, uI1),
908 gmm::sub_vector(*V, I1));
909 vars_vec_done.insert(vname1);
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),
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) {
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));
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));
933 gmm::add(M, gmm::sub_matrix(*KQJpr, I1, I2));
936 vars_mat_done.insert(std::make_pair(vname1,vname2));
945 void ga_workspace::set_include_empty_int_points(
bool include) {
946 include_empty_int_pts = include;
949 bool ga_workspace::include_empty_int_points()
const {
950 return include_empty_int_pts;
953 void ga_workspace::add_temporary_interval_for_unreduced_variable
954 (
const std::string &name)
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()) {
963 tmp_var_intervals[name] = gmm::sub_interval(nb_tmp_dof, nd);
969 void ga_workspace::clear_expressions() { trees.clear(); }
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);
981 void ga_workspace::tree_description::copy(
const tree_description& td) {
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;
993 if (td.ptree) ptree =
new ga_tree(*(td.ptree));
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; }
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())
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) {
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;
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;
1035 first_intern_dof = nb_prim_dof;
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())
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;
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)
1053 ga_workspace::~ga_workspace() { clear_expressions(); }
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)+
")";
1074 return constant_term;
1081 std::string ga_workspace::extract_order0_term() {
1083 for (
const ga_workspace::tree_description &td : trees) {
1084 if (td.order == 0) {
1085 ga_tree &local_tree = *(td.ptree);
1087 term +=
"+("+ga_tree_to_string(local_tree)+
")";
1089 term =
"("+ga_tree_to_string(local_tree)+
")";
1100 std::string ga_workspace::extract_order1_term(
const std::string &varname) {
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);
1106 term +=
"+("+ga_tree_to_string(local_tree)+
")";
1108 term =
"("+ga_tree_to_string(local_tree)+
")";
1118 std::string ga_workspace::extract_Neumann_term(
const std::string &varname) {
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);
`‘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)
*/
void copy(const L1 &l1, L2 &l2)
*/
void fill(L &l, typename gmm::linalg_traits< L >::value_type x)
*/
void clear(L &l)
clear (fill with zeros) a vector or matrix.
void resize(V &v, size_type n)
*/
void mult(const L1 &l1, const L2 &l2, L3 &l3)
*/
void add(const L1 &l1, L2 &l2)
*/
size_t size_type
used as the common size type in the library
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.