23 #include "getfem/getfem_generic_assembly_compile_and_exec.h"
24 #include "getfem/getfem_generic_assembly_functions_and_operators.h"
34 void ga_workspace::init() {
37 K = std::make_shared<model_real_sparse_matrix>(2,2);
38 V = std::make_shared<base_vector>(2);
39 KQJpr = std::make_shared<model_real_sparse_matrix>(2,2);
41 add_interpolate_transformation
43 add_interpolate_transformation
47 pstring s1 = std::make_shared<std::string>(
"Hess_u");
48 tree1.add_name(s1->c_str(), 6, 0, s1);
49 tree1.root->name =
"u";
50 tree1.root->op_type = GA_NAME;
51 tree1.root->node_type = GA_NODE_MACRO_PARAM;
53 tree1.root->nbc2 = ga_parse_prefix_operator(*s1);
54 tree1.root->nbc3 = ga_parse_prefix_test(*s1);
55 ga_macro gam1(
"Hess", tree1, 1);
56 macro_dict.add_macro(gam1);
59 pstring s2 = std::make_shared<std::string>(
"Div_u");
60 tree2.add_name(s2->c_str(), 5, 0, s2);
61 tree2.root->name =
"u";
62 tree2.root->op_type = GA_NAME;
63 tree2.root->node_type = GA_NODE_MACRO_PARAM;
65 tree2.root->nbc2 = ga_parse_prefix_operator(*s2);
66 tree2.root->nbc3 = ga_parse_prefix_test(*s2);
67 ga_macro gam2(
"Div", tree2, 1);
68 macro_dict.add_macro(gam2);
72 void ga_workspace::add_fem_variable
73 (
const std::string &name,
const mesh_fem &mf,
74 const gmm::sub_interval &I,
const model_real_plain_vector &VV) {
75 GMM_ASSERT1(nb_intern_dof == 0 || I.last() < first_intern_dof,
76 "The provided interval overlaps with internal dofs");
77 nb_prim_dof = std::max(nb_prim_dof, I.last());
78 variables.emplace(name, var_description(
true, &mf, 0, I, &VV, 1));
81 void ga_workspace::add_im_variable
82 (
const std::string &name,
const im_data &imd,
83 const gmm::sub_interval &I,
const model_real_plain_vector &VV) {
84 GMM_ASSERT1(nb_intern_dof == 0 || I.last() <= first_intern_dof,
85 "The provided interval overlaps with internal dofs");
86 nb_prim_dof = std::max(nb_prim_dof, I.last());
87 variables.emplace(name, var_description(
true, 0, &imd, I, &VV, 1));
90 void ga_workspace::add_internal_im_variable
91 (
const std::string &name,
const im_data &imd,
92 const gmm::sub_interval &I,
const model_real_plain_vector &VV) {
93 GMM_ASSERT1(I.first() >= nb_prim_dof,
94 "The provided interval overlaps with primary dofs");
95 nb_intern_dof += first_intern_dof - std::min(first_intern_dof, I.first());
96 first_intern_dof = std::min(first_intern_dof, I.first());
97 nb_intern_dof += first_intern_dof + nb_intern_dof
98 - std::min(first_intern_dof + nb_intern_dof, I.last());
99 variables.emplace(name, var_description(
true, 0, &imd, I, &VV, 1,
true));
102 void ga_workspace::add_fixed_size_variable
103 (
const std::string &name,
104 const gmm::sub_interval &I,
const model_real_plain_vector &VV) {
105 GMM_ASSERT1(nb_intern_dof == 0 || I.last() <= first_intern_dof,
106 "The provided interval overlaps with internal dofs");
107 nb_prim_dof = std::max(nb_prim_dof, I.last());
108 variables.emplace(name, var_description(
true, 0, 0, I, &VV,
109 dim_type(gmm::vect_size(VV))));
112 void ga_workspace::add_fem_constant
113 (
const std::string &name,
const mesh_fem &mf,
114 const model_real_plain_vector &VV) {
115 GMM_ASSERT1(mf.nb_dof(),
"The provided mesh_fem of variable" << name
116 <<
"has zero degrees of freedom");
117 size_type Q = gmm::vect_size(VV)/mf.nb_dof();
119 variables.emplace(name, var_description(
false, &mf, 0,
120 gmm::sub_interval(), &VV, Q));
123 void ga_workspace::add_fixed_size_constant
124 (
const std::string &name,
const model_real_plain_vector &VV) {
125 variables.emplace(name, var_description(
false, 0, 0,
126 gmm::sub_interval(), &VV,
127 gmm::vect_size(VV)));
130 void ga_workspace::add_im_data(
const std::string &name,
const im_data &imd,
131 const model_real_plain_vector &VV) {
132 variables.emplace(name, var_description
133 (
false, 0, &imd, gmm::sub_interval(), &VV,
134 gmm::vect_size(VV)/(imd.nb_filtered_index() * imd.nb_tensor_elem())));
137 bool ga_workspace::is_internal_variable(
const std::string &name)
const {
139 if ((md && md->variable_exists(name) && md->is_internal_variable(name)) ||
140 (parent_workspace && parent_workspace->variable_exists(name)
141 && parent_workspace->is_internal_variable(name)))
144 VAR_SET::const_iterator it = variables.find(name);
145 return it == variables.end() ? false : it->second.is_internal;
149 bool ga_workspace::variable_exists(
const std::string &name)
const {
150 return (md && md->variable_exists(name)) ||
151 (parent_workspace && parent_workspace->variable_exists(name)) ||
152 (variables.find(name) != variables.end());
155 bool ga_workspace::variable_group_exists(
const std::string &name)
const {
156 return (variable_groups.find(name) != variable_groups.end()) ||
157 (md && md->variable_group_exists(name)) ||
158 (parent_workspace && parent_workspace->variable_group_exists(name));
161 const std::vector<std::string>&
162 ga_workspace::variable_group(
const std::string &group_name)
const {
163 std::map<std::string, std::vector<std::string> >::const_iterator
164 it = variable_groups.find(group_name);
165 if (it != variable_groups.end())
166 return (variable_groups.find(group_name))->second;
167 if (md && md->variable_group_exists(group_name))
168 return md->variable_group(group_name);
169 if (parent_workspace &&
170 parent_workspace->variable_group_exists(group_name))
171 return parent_workspace->variable_group(group_name);
172 GMM_ASSERT1(
false,
"Undefined variable group " << group_name);
176 ga_workspace::first_variable_of_group(
const std::string &name)
const {
177 const std::vector<std::string> &t = variable_group(name);
178 GMM_ASSERT1(t.size(),
"Variable group " << name <<
" is empty");
182 bool ga_workspace::is_constant(
const std::string &name)
const {
183 VAR_SET::const_iterator it = variables.find(name);
184 if (it != variables.end())
185 return !(it->second.is_variable);
186 else if (variable_group_exists(name))
187 return is_constant(first_variable_of_group(name));
188 else if (reenabled_var_intervals.count(name))
190 else if (md && md->variable_exists(name))
191 return md->is_data(name);
192 else if (parent_workspace && parent_workspace->variable_exists(name))
193 return parent_workspace->is_constant(name);
194 GMM_ASSERT1(
false,
"Undefined variable " << name);
197 bool ga_workspace::is_disabled_variable(
const std::string &name)
const {
198 VAR_SET::const_iterator it = variables.find(name);
199 if (it != variables.end())
201 else if (variable_group_exists(name))
202 return is_disabled_variable(first_variable_of_group(name));
203 else if (reenabled_var_intervals.count(name))
205 else if (md && md->variable_exists(name))
206 return md->is_disabled_variable(name);
207 else if (parent_workspace && parent_workspace->variable_exists(name))
208 return parent_workspace->is_disabled_variable(name);
209 GMM_ASSERT1(
false,
"Undefined variable " << name);
213 ga_workspace::factor_of_variable(
const std::string &name)
const {
214 static const scalar_type one(1);
215 VAR_SET::const_iterator it = variables.find(name);
216 if (it != variables.end())
return one;
217 if (variable_group_exists(name))
219 if (md && md->variable_exists(name))
return md->factor_of_variable(name);
220 if (parent_workspace && parent_workspace->variable_exists(name))
221 return parent_workspace->factor_of_variable(name);
222 GMM_ASSERT1(
false,
"Undefined variable " << name);
225 const gmm::sub_interval &
226 ga_workspace::interval_of_variable(
const std::string &name)
const {
227 VAR_SET::const_iterator it = variables.find(name);
228 if (it != variables.end())
return it->second.I;
229 const auto it2 = reenabled_var_intervals.find(name);
230 if (it2 != reenabled_var_intervals.end())
return it2->second;
231 if (with_parent_variables && md && md->variable_exists(name))
232 return md->interval_of_variable(name);
233 else if (with_parent_variables &&
234 parent_workspace && parent_workspace->variable_exists(name))
235 return parent_workspace->interval_of_variable(name);
236 GMM_ASSERT1(
false,
"Undefined variable " << name);
240 ga_workspace::associated_mf(
const std::string &name)
const {
241 VAR_SET::const_iterator it = variables.find(name);
242 if (it != variables.end())
243 return it->second.mf;
244 if (md && md->variable_exists(name))
245 return md->pmesh_fem_of_variable(name);
246 if (parent_workspace && parent_workspace->variable_exists(name))
247 return parent_workspace->associated_mf(name);
248 if (variable_group_exists(name))
249 return associated_mf(first_variable_of_group(name));
250 GMM_ASSERT1(
false,
"Undefined variable or group " << name);
254 ga_workspace::associated_im_data(
const std::string &name)
const {
255 VAR_SET::const_iterator it = variables.find(name);
256 if (it != variables.end())
return it->second.imd;
257 if (md && md->variable_exists(name))
258 return md->pim_data_of_variable(name);
259 if (parent_workspace && parent_workspace->variable_exists(name))
260 return parent_workspace->associated_im_data(name);
261 if (variable_group_exists(name))
return 0;
262 GMM_ASSERT1(
false,
"Undefined variable " << name);
265 size_type ga_workspace::qdim(
const std::string &name)
const {
266 VAR_SET::const_iterator it = variables.find(name);
267 if (it != variables.end()) {
268 const mesh_fem *mf = it->second.mf;
269 const im_data *imd = it->second.imd;
272 return n * mf->get_qdim();
274 return n * imd->tensor_size().total_size();
278 if (md && md->variable_exists(name))
279 return md->qdim_of_variable(name);
280 if (parent_workspace && parent_workspace->variable_exists(name))
281 return parent_workspace->qdim(name);
282 if (variable_group_exists(name))
283 return qdim(first_variable_of_group(name));
284 GMM_ASSERT1(
false,
"Undefined variable or group " << name);
288 ga_workspace::qdims(
const std::string &name)
const {
289 VAR_SET::const_iterator it = variables.find(name);
290 if (it != variables.end()) {
291 const mesh_fem *mf = it->second.mf;
292 const im_data *imd = it->second.imd;
295 bgeot::multi_index mi = mf->get_qdims();
296 if (n > 1 || it->second.qdims.size() > 1) {
298 if (mi.back() == 1) { mi.back() *= it->second.qdims[0]; ++i; }
299 for (; i < it->second.qdims.size(); ++i)
300 mi.push_back(it->second.qdims[i]);
304 bgeot::multi_index mi = imd->tensor_size();
305 size_type q = n / imd->nb_filtered_index();
306 GMM_ASSERT1(q % imd->nb_tensor_elem() == 0,
307 "Invalid mesh im data vector");
308 if (n > 1 || it->second.qdims.size() > 1) {
310 if (mi.back() == 1) { mi.back() *= it->second.qdims[0]; ++i; }
311 for (; i < it->second.qdims.size(); ++i)
312 mi.push_back(it->second.qdims[i]);
316 return it->second.qdims;
318 if (md && md->variable_exists(name))
319 return md->qdims_of_variable(name);
320 if (parent_workspace && parent_workspace->variable_exists(name))
321 return parent_workspace->qdims(name);
322 if (variable_group_exists(name))
323 return qdims(first_variable_of_group(name));
324 GMM_ASSERT1(
false,
"Undefined variable or group " << name);
327 const model_real_plain_vector &
328 ga_workspace::value(
const std::string &name)
const {
329 VAR_SET::const_iterator it = variables.find(name);
330 if (it != variables.end())
331 return *(it->second.V);
332 if (md && md->variable_exists(name))
333 return md->real_variable(name);
334 if (parent_workspace && parent_workspace->variable_exists(name))
335 return parent_workspace->value(name);
336 if (variable_group_exists(name))
337 return value(first_variable_of_group(name));
338 GMM_ASSERT1(
false,
"Undefined variable or group " << name);
341 scalar_type ga_workspace::get_time_step()
const {
342 if (md)
return md->get_time_step();
343 if (parent_workspace)
return parent_workspace->get_time_step();
344 GMM_ASSERT1(
false,
"No time step defined here");
347 void ga_workspace::add_interpolate_transformation
348 (
const std::string &name, pinterpolate_transformation ptrans) {
349 if (secondary_domain_exists(name))
350 GMM_ASSERT1(
false,
"A secondary domain with the same "
351 "name already exists");
352 if (transformations.find(name) != transformations.end())
353 GMM_ASSERT1(name !=
"neighbor_element",
"neighbor_element is a "
354 "reserved interpolate transformation name");
355 transformations[name] = ptrans;
358 bool ga_workspace::interpolate_transformation_exists
359 (
const std::string &name)
const {
360 return (md && md->interpolate_transformation_exists(name)) ||
362 parent_workspace->interpolate_transformation_exists(name)) ||
363 (transformations.find(name) != transformations.end());
366 pinterpolate_transformation
367 ga_workspace::interpolate_transformation(
const std::string &name)
const {
368 auto it = transformations.find(name);
369 if (it != transformations.end())
return it->second;
370 if (md && md->interpolate_transformation_exists(name))
371 return md->interpolate_transformation(name);
372 if (parent_workspace &&
373 parent_workspace->interpolate_transformation_exists(name))
374 return parent_workspace->interpolate_transformation(name);
375 GMM_ASSERT1(
false,
"Inexistent transformation " << name);
378 bool ga_workspace::elementary_transformation_exists
379 (
const std::string &name)
const {
380 return (md && md->elementary_transformation_exists(name)) ||
382 parent_workspace->elementary_transformation_exists(name)) ||
383 (elem_transformations.find(name) != elem_transformations.end());
386 pelementary_transformation
387 ga_workspace::elementary_transformation(
const std::string &name)
const {
388 auto it = elem_transformations.find(name);
389 if (it != elem_transformations.end())
return it->second;
390 if (md && md->elementary_transformation_exists(name))
391 return md->elementary_transformation(name);
392 if (parent_workspace &&
393 parent_workspace->elementary_transformation_exists(name))
394 return parent_workspace->elementary_transformation(name);
395 GMM_ASSERT1(
false,
"Inexistent elementary transformation " << name);
398 void ga_workspace::add_secondary_domain(
const std::string &name,
399 psecondary_domain psecdom) {
400 if (interpolate_transformation_exists(name))
401 GMM_ASSERT1(
false,
"An interpolate transformation with the same "
402 "name already exists");
403 secondary_domains[name] = psecdom;
406 bool ga_workspace::secondary_domain_exists
407 (
const std::string &name)
const {
408 return (md && md->secondary_domain_exists(name)) ||
410 parent_workspace->secondary_domain_exists(name)) ||
411 (secondary_domains.find(name) != secondary_domains.end());
415 ga_workspace::secondary_domain(
const std::string &name)
const {
416 auto it = secondary_domains.find(name);
417 if (it != secondary_domains.end())
return it->second;
418 if (md && md->secondary_domain_exists(name))
419 return md->secondary_domain(name);
420 if (parent_workspace &&
421 parent_workspace->secondary_domain_exists(name))
422 return parent_workspace->secondary_domain(name);
423 GMM_ASSERT1(
false,
"Inexistent secondary domain " << name);
428 ga_workspace::register_region(
const mesh &m,
const mesh_region ®ion) {
429 if (&m == &dummy_mesh())
432 std::list<mesh_region> &lmr = registered_mesh_regions[&m];
433 for (
const mesh_region &rg : lmr)
434 if (rg.compare(m, region, m))
return rg;
435 lmr.push_back(region);
439 void ga_workspace::add_tree(ga_tree &tree,
const mesh &m,
440 const mesh_im &mim,
const mesh_region &rg,
442 bool function_expr, operation_type op_type,
443 const std::string varname_interpolation) {
450 if ((tree.root->test_function_type >= 1 &&
451 is_disabled_variable(tree.root->name_test1)) ||
452 (tree.root->test_function_type >= 2 &&
453 is_disabled_variable(tree.root->name_test2))) {
461 if (op_type != ga_workspace::ASSEMBLY)
462 order = add_derivative_order;
464 switch(tree.root->test_function_type) {
465 case 0: order = 0;
break;
466 case 1: order = 1;
break;
467 case 3: order = 2;
break;
468 default: GMM_ASSERT1(
false,
"Inconsistent term "
469 << tree.root->test_function_type);
474 for (
const ga_workspace::tree_description &td : trees) {
475 if (td.mim == &mim &&
477 td.secondary_domain == tree.secondary_domain &&
479 td.name_test1 == tree.root->name_test1 &&
480 td.interpolate_name_test1 == tree.root->interpolate_name_test1 &&
481 td.name_test2 == tree.root->name_test2 &&
482 td.interpolate_name_test2 == tree.root->interpolate_name_test2 &&
484 td.operation == op_type &&
485 td.varname_interpolation == varname_interpolation) {
486 ga_tree &ftree = *(td.ptree);
488 ftree.insert_node(ftree.root, GA_NODE_OP);
489 ftree.root->op_type = GA_PLUS;
490 ftree.root->children.resize(2,
nullptr);
491 ftree.copy_node(tree.root, ftree.root->children[1]);
492 ftree.root->children[1]->parent = ftree.root;
493 ga_semantic_analysis(ftree, *
this, m,
494 ref_elt_dim_of_mesh(m,rg),
false, function_expr);
501 ind_tree = trees.size();
503 trees.push_back(tree_description());
504 trees.back().mim = &mim;
506 trees.back().rg = &rg;
507 trees.back().secondary_domain = tree.secondary_domain;
508 trees.back().ptree =
new ga_tree;
509 trees.back().ptree->swap(tree);
510 pga_tree_node root = trees.back().ptree->root;
511 trees.back().name_test1 = root->name_test1;
512 trees.back().name_test2 = root->name_test2;
513 trees.back().interpolate_name_test1 = root->interpolate_name_test1;
514 trees.back().interpolate_name_test2 = root->interpolate_name_test2;
515 trees.back().order = order;
516 trees.back().operation = op_type;
517 trees.back().varname_interpolation = varname_interpolation;
520 if (op_type == ga_workspace::ASSEMBLY && order < add_derivative_order) {
521 std::set<var_trans_pair> expr_variables;
522 ga_extract_variables((remain ? tree : *(trees[ind_tree].ptree)).root,
523 *
this, m, expr_variables,
true);
524 for (
const var_trans_pair &var : expr_variables) {
525 if (!(is_constant(var.varname))) {
526 ga_tree dtree = (remain ? tree : *(trees[ind_tree].ptree));
530 ga_derivative(dtree, *
this, m, var.varname, var.transname, 1+order);
533 ga_semantic_analysis(dtree, *
this, m,
534 ref_elt_dim_of_mesh(m,rg),
false,function_expr);
537 add_tree(dtree, m, mim, rg, add_derivative_order,
538 function_expr, op_type, varname_interpolation);
545 size_type ga_workspace::add_expression(
const std::string &expr,
547 const mesh_region &rg_,
549 const std::string &secondary_dom) {
550 const mesh_region &rg = register_region(mim.linked_mesh(), rg_);
554 std::vector<ga_tree> ltrees(1);
555 ga_read_string(expr, ltrees[0], macro_dictionary());
556 if (secondary_dom.size()) {
557 GMM_ASSERT1(secondary_domain_exists(secondary_dom),
558 "Unknown secondary domain " << secondary_dom);
559 ltrees[0].secondary_domain = secondary_dom;
562 ga_semantic_analysis(ltrees[0], *
this, mim.linked_mesh(),
563 ref_elt_dim_of_mesh(mim.linked_mesh(),rg),
566 GA_TOC(
"First analysis time");
567 if (ltrees[0].root) {
568 if (test1.size() > 1 || test2.size() > 1) {
571 test2.insert(var_trans_pair());
572 ltrees.resize(test1.size()*test2.size(), ltrees[0]);
573 auto ltree = ltrees.begin();
574 for (
const auto &t1 : test1) {
575 for (
const auto &t2 : test2) {
577 if (ntest2 > 0) selected_test2 = t2;
579 ga_semantic_analysis(*ltree, *
this, mim.linked_mesh(),
580 ref_elt_dim_of_mesh(mim.linked_mesh(),rg),
583 if (ltree != ltrees.end()) ++ltree;
586 if (ntest2 == 0) test2.clear();
589 for (ga_tree <ree : ltrees) {
592 max_order = std::max(ltree.root->nb_test_functions(), max_order);
593 add_tree(ltree, mim.linked_mesh(), mim, rg,
594 add_derivative_order,
true);
598 GA_TOC(
"Time for add expression");
602 void ga_workspace::add_function_expression(
const std::string &expr) {
604 ga_read_string(expr, tree, macro_dictionary());
605 ga_semantic_analysis(tree, *
this, dummy_mesh(), 1,
false,
true);
614 void ga_workspace::add_interpolation_expression(
const std::string &expr,
616 const mesh_region &rg_) {
617 const mesh_region &rg = register_region(m, rg_);
619 ga_read_string(expr, tree, macro_dictionary());
620 ga_semantic_analysis(tree, *
this, m, ref_elt_dim_of_mesh(m,rg),
626 ga_workspace::PRE_ASSIGNMENT);
630 void ga_workspace::add_interpolation_expression(
const std::string &expr,
632 const mesh_region &rg_) {
633 const mesh &m = mim.linked_mesh();
634 const mesh_region &rg = register_region(m, rg_);
636 ga_read_string(expr, tree, macro_dictionary());
637 ga_semantic_analysis(tree, *
this, m, ref_elt_dim_of_mesh(m,rg),
640 GMM_ASSERT1(tree.root->nb_test_functions() == 0,
641 "Invalid expression containing test functions");
642 add_tree(tree, m, mim, rg, 0,
false,
643 ga_workspace::PRE_ASSIGNMENT);
647 void ga_workspace::add_assignment_expression
648 (
const std::string &varname,
const std::string &expr,
const mesh_region &rg_,
650 const im_data *imd = associated_im_data(varname);
651 GMM_ASSERT1(imd != 0,
"Only applicable to im_data");
652 const mesh_im &mim = imd->linked_mesh_im();
653 const mesh &m = mim.linked_mesh();
654 const mesh_region &rg = register_region(m, rg_);
656 ga_read_string(expr, tree, macro_dictionary());
657 ga_semantic_analysis(tree, *
this, m, ref_elt_dim_of_mesh(m,rg),
false,
false);
659 GMM_ASSERT1(tree.root->nb_test_functions() == 0,
660 "Invalid expression containing test functions");
661 add_tree(tree, m, mim, rg, order,
false,
662 before ? ga_workspace::PRE_ASSIGNMENT
663 : ga_workspace::POST_ASSIGNMENT,
668 bool ga_workspace::used_variables(std::vector<std::string> &vl,
669 std::vector<std::string> &vl_test1,
670 std::vector<std::string> &vl_test2,
671 std::vector<std::string> &dl,
674 std::set<var_trans_pair> vll, dll;
675 for (
const std::string &v : vl) vll.insert(var_trans_pair(v,
""));
676 for (
const std::string &d : dl) dll.insert(var_trans_pair(d,
""));
678 for (
const ga_workspace::tree_description &td : trees) {
679 std::set<var_trans_pair> dllaux;
680 bool fv = ga_extract_variables(td.ptree->root, *
this, *(td.m),
683 if (td.order == order)
684 for (
const auto &t : dllaux)
690 if (td.order == order) {
691 if (variable_group_exists(td.name_test1)) {
692 for (
const std::string &t : variable_group(td.name_test1))
693 vll.insert(var_trans_pair(t, td.interpolate_name_test1));
695 vll.insert(var_trans_pair(td.name_test1,
696 td.interpolate_name_test1));
699 for (
const std::string &t : vl_test1)
700 if (td.name_test1 == t)
703 vl_test1.push_back(td.name_test1);
707 if (td.order == order) {
708 if (variable_group_exists(td.name_test1)) {
709 for (
const std::string &t : variable_group(td.name_test1))
710 vll.insert(var_trans_pair(t, td.interpolate_name_test1));
712 vll.insert(var_trans_pair(td.name_test1,
713 td.interpolate_name_test1));
715 if (variable_group_exists(td.name_test2)) {
716 for (
const std::string &t : variable_group(td.name_test2))
717 vll.insert(var_trans_pair(t, td.interpolate_name_test2));
719 vll.insert(var_trans_pair(td.name_test2,
720 td.interpolate_name_test2));
723 for (
size_type j = 0; j < vl_test1.size(); ++j)
724 if ((td.name_test1 == vl_test1[j]) &&
725 (td.name_test2 == vl_test2[j]))
728 vl_test1.push_back(td.name_test1);
729 vl_test2.push_back(td.name_test2);
732 if (fv) islin =
false;
737 for (
const auto &var : vll)
738 if (vl.size() == 0 || var.varname != vl.back())
739 vl.push_back(var.varname);
741 for (
const auto &var : dll)
742 if (dl.size() == 0 || var.varname != dl.back())
743 dl.push_back(var.varname);
748 bool ga_workspace::is_linear(
size_type order) {
749 std::vector<std::string> vl, vl_test1, vl_test2, dl;
750 return used_variables(vl, vl_test1, vl_test2, dl, order);
753 void ga_workspace::define_variable_group(
const std::string &group_name,
754 const std::vector<std::string> &nl) {
755 GMM_ASSERT1(!(variable_exists(group_name)),
"The name of a group of "
756 "variables cannot be the same as a variable name");
758 std::set<const mesh *> ms;
759 bool is_data_ =
false;
760 for (
size_type i = 0; i < nl.size(); ++i) {
762 is_data_ = is_constant(nl[i]);
764 GMM_ASSERT1(is_data_ == is_constant(nl[i]),
765 "It is not possible to mix variables and data in a group");
767 GMM_ASSERT1(variable_exists(nl[i]),
768 "All variables in a group have to exist in the model");
769 const mesh_fem *mf = associated_mf(nl[i]);
770 GMM_ASSERT1(mf,
"Variables in a group should be fem variables");
771 GMM_ASSERT1(ms.find(&(mf->linked_mesh())) == ms.end(),
772 "Two variables in a group cannot share the same mesh");
773 ms.insert(&(mf->linked_mesh()));
775 variable_groups[group_name] = nl;
779 const std::string &ga_workspace::variable_in_group
780 (
const std::string &group_name,
const mesh &m)
const {
781 if (variable_group_exists(group_name)) {
782 for (
const std::string &t : variable_group(group_name))
783 if (&(associated_mf(t)->linked_mesh()) == &m)
785 GMM_ASSERT1(
false,
"No variable in this group for the given mesh");
791 void ga_workspace::assembly(
size_type order,
bool condensation) {
793 const ga_workspace *w =
this;
794 while (w->parent_workspace) w = w->parent_workspace;
795 if (w->md) w->md->nb_dof();
798 ga_instruction_set gis;
799 ga_compile(*
this, gis, order, condensation);
800 GA_TOCTIC(
"Compile time");
802 size_type nb_tot_dof = condensation ? nb_prim_dof + nb_intern_dof
813 if (KQJpr.use_count()) {
817 }
else if (condensation)
818 GMM_ASSERT1(gmm::mat_nrows(*KQJpr) == nb_tot_dof &&
819 gmm::mat_ncols(*KQJpr) == nb_prim_dof,
"Wrong sizes");
823 gmm::resize(col_unreduced_K, nb_tot_dof, nb_tmp_dof);
824 gmm::resize(row_unreduced_K, nb_tmp_dof, nb_prim_dof);
825 gmm::resize(row_col_unreduced_K, nb_tmp_dof, nb_tmp_dof);
832 if (order == 1 || (order == 2 && condensation)) {
833 if (order == 2 && condensation) {
834 GMM_ASSERT1(V->size() == nb_tot_dof,
835 "Wrong size of assembled vector in workspace");
839 }
else if (V.use_count()) {
843 GMM_ASSERT1(V->size() == nb_tot_dof,
844 "Wrong size of assembled vector in workspace");
850 GA_TOCTIC(
"Init time");
852 GA_TOCTIC(
"Exec time");
855 MPI_SUM_VECTOR(assemb_t.as_vector());
856 }
else if (order == 1 || (order == 2 && condensation)) {
858 MPI_SUM_VECTOR(unreduced_V);
864 std::set<std::string> vars_vec_done;
865 std::set<std::pair<std::string, std::string> > vars_mat_done;
866 for (
const auto &term : gis.unreduced_terms) {
867 const std::string &name1 = term.first;
868 const std::string &name2 = term.second;
869 const std::vector<std::string>
870 vg1_(1,name1), vg2_(1,name2),
871 &vg1 = variable_group_exists(name1) ? variable_group(name1) : vg1_,
872 &vg2 = variable_group_exists(name2) ? variable_group(name2) : vg2_;
874 for (
const std::string &vname1 : vg1) {
875 const mesh_fem *mf1 = associated_mf(vname1);
876 if (mf1 && mf1->is_reduced() && vars_vec_done.count(vname1) == 0) {
877 gmm::sub_interval uI1 = temporary_interval_of_variable(vname1),
878 I1 = interval_of_variable(vname1);
880 gmm::sub_vector(unreduced_V, uI1),
881 gmm::sub_vector(*V, I1));
882 vars_vec_done.insert(vname1);
886 for (
const std::string &vname1 : vg1) {
887 for (
const std::string &vname2 : vg2) {
888 const mesh_fem *mf1 = associated_mf(vname1),
889 *mf2 = associated_mf(vname2);
890 if (((mf1 && mf1->is_reduced()) || (mf2 && mf2->is_reduced()))
891 && vars_mat_done.count(std::make_pair(vname1,vname2)) == 0) {
893 uI1 = temporary_interval_of_variable(vname1),
894 uI2 = temporary_interval_of_variable(vname2),
895 I1 = interval_of_variable(vname1),
896 I2 = interval_of_variable(vname2);
897 if (mf1 && mf1->is_reduced()) {
898 if (condensation && vars_vec_done.count(vname1) == 0) {
900 gmm::sub_vector(unreduced_V, uI1),
901 gmm::sub_vector(*V, I1));
902 vars_vec_done.insert(vname1);
904 if (mf2 && mf2->is_reduced()) {
905 model_real_sparse_matrix aux(I1.size(), uI2.size());
906 model_real_row_sparse_matrix M(I1.size(), I2.size());
907 gmm::mult(gmm::transposed(mf1->extension_matrix()),
908 gmm::sub_matrix(row_col_unreduced_K, uI1, uI2),
910 gmm::mult(aux, mf2->extension_matrix(), M);
911 gmm::add(M, gmm::sub_matrix(*K, I1, I2));
912 }
else if (I2.first() < nb_prim_dof) {
913 model_real_sparse_matrix M(I1.size(), I2.size());
914 gmm::mult(gmm::transposed(mf1->extension_matrix()),
915 gmm::sub_matrix(row_unreduced_K, uI1, I2), M);
916 gmm::add(M, gmm::sub_matrix(*K, I1, I2));
919 model_real_row_sparse_matrix M(I1.size(), I2.size());
920 gmm::mult(gmm::sub_matrix(col_unreduced_K, I1, uI2),
921 mf2->extension_matrix(), M);
922 if (I1.first() < nb_prim_dof) {
923 GMM_ASSERT1(I1.last() <= nb_prim_dof,
"Internal error");
924 gmm::add(M, gmm::sub_matrix(*K, I1, I2));
926 gmm::add(M, gmm::sub_matrix(*KQJpr, I1, I2));
929 vars_mat_done.insert(std::make_pair(vname1,vname2));
938 void ga_workspace::set_include_empty_int_points(
bool include) {
939 include_empty_int_pts = include;
942 bool ga_workspace::include_empty_int_points()
const {
943 return include_empty_int_pts;
946 void ga_workspace::add_temporary_interval_for_unreduced_variable
947 (
const std::string &name)
949 if (variable_group_exists(name)) {
950 for (
const std::string &v : variable_group(name))
951 add_temporary_interval_for_unreduced_variable(v);
952 }
else if (tmp_var_intervals.count(name) == 0) {
953 const mesh_fem *mf = associated_mf(name);
954 if (mf && mf->is_reduced()) {
956 tmp_var_intervals[name] = gmm::sub_interval(nb_tmp_dof, nd);
962 void ga_workspace::clear_expressions() { trees.clear(); }
964 void ga_workspace::print(std::ostream &str) {
965 for (
size_type i = 0; i < trees.size(); ++i)
966 if (trees[i].ptree->root) {
967 cout <<
"Expression tree " << i <<
" of order " <<
968 trees[i].ptree->root->nb_test_functions() <<
" :" << endl;
969 ga_print_node(trees[i].ptree->root, str);
974 void ga_workspace::tree_description::copy(
const tree_description& td) {
976 operation = td.operation;
977 varname_interpolation = td.varname_interpolation;
978 name_test1 = td.name_test1;
979 name_test2 = td.name_test2;
980 interpolate_name_test1 = td.interpolate_name_test1;
981 interpolate_name_test2 = td.interpolate_name_test2;
986 if (td.ptree) ptree =
new ga_tree(*(td.ptree));
989 ga_workspace::tree_description &ga_workspace::tree_description::operator =
990 (
const ga_workspace::tree_description& td)
991 {
if (ptree)
delete ptree; ptree = 0;
copy(td);
return *
this; }
992 ga_workspace::tree_description::~tree_description()
993 {
if (ptree)
delete ptree; ptree = 0; }
996 const inherit var_inherit)
997 : md(&md_), parent_workspace(0),
998 with_parent_variables(var_inherit == inherit::ENABLED ||
999 var_inherit == inherit::ALL),
1000 nb_tmp_dof(0), macro_dict(md_.macro_dictionary())
1003 nb_prim_dof = with_parent_variables ? md->nb_primary_dof() : 0;
1004 nb_intern_dof = with_parent_variables ? md->nb_internal_dof() : 0;
1005 if (var_inherit == inherit::ALL) {
1006 model::varnamelist vlmd;
1007 md->variable_list(vlmd);
1008 for (
const auto &varname : vlmd)
1009 if (md->is_disabled_variable(varname)) {
1010 if (md->is_affine_dependent_variable(varname)) {
1011 std::string orgvarname = md->org_variable(varname);
1012 if (reenabled_var_intervals.count(orgvarname) == 0) {
1013 size_type varsize = gmm::vect_size(md->real_variable(orgvarname));
1014 reenabled_var_intervals[orgvarname]
1015 = gmm::sub_interval (nb_prim_dof, varsize);
1016 nb_prim_dof += varsize;
1018 reenabled_var_intervals[varname]
1019 = reenabled_var_intervals[orgvarname];
1020 }
else if (reenabled_var_intervals.count(varname) == 0) {
1021 size_type varsize = gmm::vect_size(md->real_variable(varname));
1022 reenabled_var_intervals[varname]
1023 = gmm::sub_interval(nb_prim_dof, varsize);
1024 nb_prim_dof += varsize;
1028 first_intern_dof = nb_prim_dof;
1030 ga_workspace::ga_workspace(
const ga_workspace &gaw,
1031 const inherit var_inherit)
1032 : md(0), parent_workspace(&gaw),
1033 with_parent_variables(var_inherit == inherit::ENABLED ||
1034 var_inherit == inherit::ALL),
1035 nb_tmp_dof(0), macro_dict(gaw.macro_dictionary())
1038 nb_prim_dof = with_parent_variables ? gaw.nb_primary_dof() : 0;
1039 nb_intern_dof = with_parent_variables ? gaw.nb_internal_dof() : 0;
1040 first_intern_dof = with_parent_variables ? gaw.first_internal_dof() : 0;
1042 ga_workspace::ga_workspace()
1043 : md(0), parent_workspace(0), with_parent_variables(false),
1044 nb_prim_dof(0), nb_intern_dof(0), first_intern_dof(0), nb_tmp_dof(0)
1046 ga_workspace::~ga_workspace() { clear_expressions(); }
1052 std::string ga_workspace::extract_constant_term(
const mesh &m) {
1053 std::string constant_term;
1054 for (
const ga_workspace::tree_description &td : trees) {
1055 if (td.order == 1) {
1056 ga_tree local_tree = *(td.ptree);
1057 if (local_tree.root)
1058 ga_node_extract_constant_term(local_tree, local_tree.root, *
this, m);
1059 if (local_tree.root)
1060 ga_semantic_analysis(local_tree, *
this, m,
1061 ref_elt_dim_of_mesh(m,-1),
false,
false);
1062 if (local_tree.root && local_tree.root->node_type != GA_NODE_ZERO) {
1063 constant_term +=
"-("+ga_tree_to_string(local_tree)+
")";
1067 return constant_term;
1074 std::string ga_workspace::extract_order0_term() {
1076 for (
const ga_workspace::tree_description &td : trees) {
1077 if (td.order == 0) {
1078 const ga_tree &local_tree = *(td.ptree);
1080 term +=
"+("+ga_tree_to_string(local_tree)+
")";
1082 term =
"("+ga_tree_to_string(local_tree)+
")";
1093 std::string ga_workspace::extract_order1_term(
const std::string &varname) {
1095 for (
const ga_workspace::tree_description &td : trees) {
1096 if (td.order == 1 && td.name_test1 == varname) {
1097 const ga_tree &local_tree = *(td.ptree);
1099 term +=
"+("+ga_tree_to_string(local_tree)+
")";
1101 term =
"("+ga_tree_to_string(local_tree)+
")";
1111 std::string ga_workspace::extract_Neumann_term(
const std::string &varname) {
1113 for (
const ga_workspace::tree_description &td : trees) {
1114 if (td.order == 1 && td.name_test1 == varname) {
1115 const ga_tree &local_tree = *(td.ptree);
1116 if (local_tree.root)
1117 ga_extract_Neumann_term(local_tree, varname, *
this,
1118 local_tree.root, result);
`‘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.