37 model::model(
bool comp_version) {
38 init(); complex_version = comp_version;
39 is_linear_ = is_symmetric_ = is_coercive_ =
true;
41 time_integration = 0; init_step =
false; time_step = scalar_type(1);
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);
72 void model::var_description::set_size() {
74 v_num_var_iter.resize(n_iter);
75 v_num_iter.resize(n_iter);
76 size_type s = is_fem_dofs ? passociated_mf()->nb_dof()
77 : (imd ? imd->nb_filtered_index()
78 * imd->nb_tensor_elem()
83 complex_value[i].resize(s);
85 real_value[i].resize(s);
86 if (is_affine_dependent) {
88 affine_complex_value.resize(s);
90 affine_real_value.resize(s);
94 size_type model::var_description::add_temporary(gmm::uint64_type id_num) {
96 for (; nit < n_iter + n_temp_iter ; ++nit)
97 if (v_num_var_iter[nit] == id_num)
break;
98 if (nit >= n_iter + n_temp_iter) {
100 v_num_var_iter.resize(nit+1);
101 v_num_var_iter[nit] = id_num;
102 v_num_iter.resize(nit+1);
106 complex_value.resize(n_iter + n_temp_iter);
107 complex_value[nit].resize(s);
110 real_value.resize(n_iter + n_temp_iter);
111 real_value[nit].resize(s);
117 void model::var_description::clear_temporaries() {
121 complex_value.resize(n_iter);
123 real_value.resize(n_iter);
126 bool model::check_name_validity(
const std::string &name,
bool assert)
const {
128 if (variables.count(name) != 0) {
129 GMM_ASSERT1(!assert,
"Variable " << name <<
" already exists");
131 }
else if (variable_groups.count(name) != 0) {
133 name <<
" corresponds to an already existing group of "
138 name <<
" corresponds to an already existing macro");
140 }
else if (name.compare(
"X") == 0) {
141 GMM_ASSERT1(!assert,
"X is a reserved keyword of the generic "
142 "assembly language");
146 int ga_valid = ga_check_name_validity(name);
148 GMM_ASSERT1(!assert,
"Invalid variable name, corresponds to an "
149 "operator or function name of the generic assembly language");
151 }
else if (ga_valid == 2) {
152 GMM_ASSERT1(!assert,
"Invalid variable name having a reserved "
153 "prefix used by the generic assembly language");
155 }
else if (ga_valid == 3) {
156 std::string org_name = sup_previous_and_dot_to_varname(name);
157 if (org_name.size() < name.size() &&
158 variables.find(org_name) != variables.end()) {
160 "Dot and Previous are reserved prefix used for time "
161 "integration schemes");
166 bool valid = !name.empty() && isalpha(name[0]);
168 for (
size_type i = 1; i < name.size(); ++i)
169 if (!(isalnum(name[i]) || name[i] ==
'_')) valid =
false;
170 GMM_ASSERT1(!assert || valid,
171 "Illegal variable name : \"" << name <<
"\"");
176 std::string res_name = name;
177 bool valid = check_name_validity(res_name,
false);
178 for (
size_type i = 2; !valid && i < 50; ++i) {
180 m << name <<
'_' << i;
182 valid = check_name_validity(res_name,
false);
184 for (
size_type i = 2; !valid && i < 1000; ++i) {
186 m <<
"new_" << name <<
'_' << i;
188 valid = check_name_validity(res_name,
false);
190 GMM_ASSERT1(valid,
"Illegal variable name: " << name);
195 model::VAR_SET::const_iterator
196 model::find_variable(
const std::string &name)
const {
197 VAR_SET::const_iterator it = variables.find(name);
198 GMM_ASSERT1(it != variables.end(),
"Undefined variable " << name);
202 const model::var_description &
203 model::variable_description(
const std::string &name)
const {
204 return find_variable(name)->second;
207 std::string sup_previous_and_dot_to_varname(std::string v) {
208 if (!(v.compare(0, 8,
"Previous")) && (v[8] ==
'_' || v[9] ==
'_')) {
209 v = v.substr((v[8] ==
'_') ? 9 : 10);
211 if (!(v.compare(0, 3,
"Dot")) && (v[3] ==
'_' || v[4] ==
'_')) {
212 v = v.substr((v[3] ==
'_') ? 4 : 5);
219 return name.substr(0, PREFIX_OLD_LENGTH) ==
PREFIX_OLD;
223 return is_old(name) ? name.substr(PREFIX_OLD_LENGTH) : name;
227 if (
is_old(name))
return false;
228 VAR_SET::const_iterator it = find_variable(name);
229 if (!(it->second.is_variable))
return false;
230 if (it->second.is_affine_dependent)
231 it = variables.find(it->second.org_name);
232 return it->second.is_disabled;
236 if (
is_old(name))
return true;
237 VAR_SET::const_iterator it = find_variable(name);
238 if (it->second.is_affine_dependent)
239 it = variables.find(it->second.org_name);
240 return !(it->second.is_variable) || it->second.is_disabled;
244 return is_old(name) || !(variable_description(name).is_variable);
248 if (
is_old(name))
return false;
249 const auto &var_descr = variable_description(name);
250 return var_descr.is_internal && var_descr.is_enabled();
253 bool model::is_affine_dependent_variable(
const std::string &name)
const {
254 return !(
is_old(name)) && variable_description(name).is_affine_dependent;
258 model::org_variable(
const std::string &name)
const {
259 GMM_ASSERT1(is_affine_dependent_variable(name),
260 "For affine dependent variables only");
261 return variable_description(name).org_name;
265 model::factor_of_variable(
const std::string &name)
const {
266 return variable_description(name).alpha;
269 void model::set_factor_of_variable(
const std::string &name, scalar_type a) {
270 VAR_SET::iterator it = variables.find(name);
271 GMM_ASSERT1(it != variables.end(),
"Undefined variable " << name);
272 if (it->second.alpha != a) {
273 it->second.alpha = a;
274 for (
auto &v_num : it->second.v_num_data) v_num = act_counter();
278 bool model::is_im_data(
const std::string &name)
const {
283 model::pim_data_of_variable(
const std::string &name)
const {
287 const gmm::uint64_type &
288 model::version_number_of_data_variable(
const std::string &name,
290 VAR_SET::const_iterator it = find_variable(name);
291 if (niter ==
size_type(-1)) niter = it->second.default_iter;
292 return it->second.v_num_data[niter];
297 if (act_size_to_be_done) actualize_sizes();
299 return gmm::vect_size(crhs);
300 else if (with_internal && gmm::vect_size(full_rrhs))
301 return gmm::vect_size(full_rrhs);
303 return gmm::vect_size(rrhs);
306 void model::resize_global_system()
const {
309 for (
auto &&v : variables)
310 if (v.second.is_variable) {
311 if (v.second.is_disabled)
312 v.second.I = gmm::sub_interval(0,0);
313 else if (!v.second.is_affine_dependent && !v.second.is_internal) {
314 v.second.I = gmm::sub_interval(full_size, v.second.size());
315 full_size += v.second.size();
320 for (
auto &&v : variables)
321 if (v.second.is_internal && v.second.is_enabled()) {
322 v.second.I = gmm::sub_interval(full_size, v.second.size());
323 full_size += v.second.size();
326 for (
auto &&v : variables)
327 if (v.second.is_affine_dependent) {
328 v.second.I = variables.find(v.second.org_name)->second.I;
332 if (complex_version) {
341 if (full_size > primary_size) {
343 gmm::resize(internal_rTM, full_size-primary_size, primary_size);
352 for (dal::bv_visitor ib(valid_bricks); !ib.finished(); ++ib)
353 for (
const term_description &term : bricks[ib].tlist)
354 if (term.is_global) {
355 bricks[ib].terms_to_be_computed =
true;
360 void model::actualize_sizes()
const {
362 bool actualized =
false;
363 getfem::local_guard lock = locks_.get_lock();
364 if (actualized)
return;
366 act_size_to_be_done =
false;
368 std::map<std::string, std::vector<std::string> > multipliers;
369 std::set<std::string> tobedone;
385 for (
auto &&v : variables) {
386 const std::string &vname = v.first;
387 var_description &vdescr = v.second;
388 if (vdescr.is_fem_dofs && !vdescr.is_affine_dependent) {
389 if ((vdescr.filter & VDESCRFILTER_CTERM)
390 || (vdescr.filter & VDESCRFILTER_INFSUP)) {
391 VAR_SET::iterator vfilt = variables.find(vdescr.filter_var);
392 GMM_ASSERT1(vfilt != variables.end(),
"The primal variable of the"
393 " multiplier does not exist : " << vdescr.filter_var);
394 GMM_ASSERT1(vfilt->second.is_fem_dofs,
"The primal variable of "
395 "the multiplier is not a fem variable");
396 multipliers[vdescr.filter_var].push_back(vname);
397 if (vdescr.v_num < vdescr.mf->version_number() ||
398 vdescr.v_num < vfilt->second.mf->version_number()) {
399 tobedone.insert(vdescr.filter_var);
402 switch (vdescr.filter) {
403 case VDESCRFILTER_NO:
404 if (vdescr.v_num < vdescr.mf->version_number()) {
406 vdescr.v_num = act_counter();
409 case VDESCRFILTER_REGION:
410 if (vdescr.v_num < vdescr.mf->version_number()) {
412 dor = vdescr.mf->dof_on_region(vdescr.filter_region);
413 vdescr.partial_mf->adapt(dor);
415 vdescr.v_num = act_counter();
423 && vdescr.v_num < vdescr.imd->version_number()) {
425 vdescr.v_num = act_counter();
429 for (
auto &&v : variables) {
430 var_description &vdescr = v.second;
431 if (vdescr.is_fem_dofs && !(vdescr.is_affine_dependent) &&
432 ((vdescr.filter & VDESCRFILTER_CTERM)
433 || (vdescr.filter & VDESCRFILTER_INFSUP))) {
434 if (tobedone.count(vdescr.filter_var)) {
439 dal::bit_vector alldof; alldof.add(0, vdescr.mf->nb_dof());
440 vdescr.partial_mf->adapt(alldof);
442 vdescr.v_num = act_counter();
447 resize_global_system();
449 for (
const std::string &vname : tobedone) {
456 const std::vector<std::string> &mults = multipliers[vname];
457 const var_description &vdescr = variable_description(vname);
459 gmm::col_matrix< gmm::rsvector<scalar_type> > MGLOB;
460 if (mults.size() > 1) {
462 for (
const std::string &mult : mults)
463 s += variable_description(mult).mf->nb_dof();
467 std::set<size_type> glob_columns;
468 for (
const std::string &multname : mults) {
469 var_description &multdescr = variables.find(multname)->second;
475 gmm::col_matrix< gmm::rsvector<scalar_type> >
476 MM(vdescr.associated_mf().nb_dof(), multdescr.mf->nb_dof());
477 bool termadded =
false;
479 if (multdescr.filter & VDESCRFILTER_CTERM) {
481 for (dal::bv_visitor ib(valid_bricks); !ib.finished(); ++ib) {
482 const brick_description &brick = bricks[ib];
484 bool cplx =
is_complex() && brick.pbr->is_complex();
486 if (brick.tlist.size() == 0) {
487 bool varc =
false, multc =
false;
488 for (
const std::string &var : brick.vlist) {
489 if (multname.compare(var) == 0) multc =
true;
490 if (vname.compare(var) == 0) varc =
true;
493 GMM_ASSERT1(!cplx,
"Sorry, not taken into account");
494 generic_expressions.clear();
495 brick.terms_to_be_computed =
true;
496 update_brick(ib, BUILD_MATRIX);
497 if (generic_expressions.size()) {
498 GMM_TRACE2(
"Generic assembly for actualize sizes");
501 accumulated_distro<decltype(rTM)> distro_rTM(rTM);
503 ga_workspace workspace(*
this);
504 for (
const auto &ge : generic_expressions)
505 workspace.add_expression(ge.expr, ge.mim, ge.region,
506 2, ge.secondary_domain);
507 workspace.set_assembled_matrix(distro_rTM);
508 workspace.assembly(2);
511 gmm::add(gmm::sub_matrix(rTM, vdescr.I, multdescr.I), MM);
512 gmm::add(gmm::transposed
513 (gmm::sub_matrix(rTM, multdescr.I, vdescr.I)), MM);
519 for (
size_type j = 0; j < brick.tlist.size(); ++j) {
520 const term_description &term = brick.tlist[j];
522 if (term.is_matrix_term) {
523 if (term.is_global) {
524 bool varc =
false, multc =
false;
525 for (
const std::string &var : brick.vlist) {
526 if (multname.compare(var) == 0) multc =
true;
527 if (vname.compare(var) == 0) varc =
true;
530 GMM_ASSERT1(!cplx,
"Sorry, not taken into account");
531 generic_expressions.clear();
533 brick.terms_to_be_computed =
true;
534 update_brick(ib, BUILD_MATRIX);
537 gmm::add(gmm::sub_matrix(brick.rmatlist[j],
538 multdescr.I, vdescr.I),
540 gmm::add(gmm::transposed(gmm::sub_matrix
542 vdescr.I, multdescr.I)),
546 }
else if (multname.compare(term.var1) == 0 &&
547 vname.compare(term.var2) == 0) {
549 brick.terms_to_be_computed =
true;
550 update_brick(ib, BUILD_MATRIX);
554 gmm::add(gmm::transposed(gmm::real_part(brick.cmatlist[j])),
557 gmm::add(gmm::transposed(brick.rmatlist[j]), MM);
560 }
else if (multname.compare(term.var2) == 0 &&
561 vname.compare(term.var1) == 0) {
563 brick.terms_to_be_computed =
true;
564 update_brick(ib, BUILD_MATRIX);
568 gmm::add(gmm::real_part(brick.cmatlist[j]), MM);
570 gmm::add(brick.rmatlist[j], MM);
578 GMM_WARNING1(
"No term found to filter multiplier " << multname
579 <<
". Variable is cancelled");
580 }
else if (multdescr.filter & VDESCRFILTER_INFSUP) {
581 mesh_region rg(multdescr.filter_region);
582 multdescr.filter_mim->linked_mesh().intersect_with_mpi_region(rg);
584 *(multdescr.mf), rg);
587 MPI_SUM_SPARSE_MATRIX(MM);
592 std::set<size_type> columns;
594 if (columns.size() == 0)
595 GMM_WARNING1(
"Empty basis found for multiplier " << multname);
597 if (mults.size() > 1) {
598 gmm::copy(MM, gmm::sub_matrix
600 gmm::sub_interval(0, vdescr.associated_mf().nb_dof()),
601 gmm::sub_interval(s, multdescr.mf->nb_dof())));
603 glob_columns.insert(s + icol);
604 s += multdescr.mf->nb_dof();
606 dal::bit_vector kept;
609 if (multdescr.filter & VDESCRFILTER_REGION)
610 kept &= multdescr.mf->dof_on_region(multdescr.filter_region);
611 multdescr.partial_mf->adapt(kept);
612 multdescr.set_size();
613 multdescr.v_num = act_counter();
621 if (mults.size() > 1) {
629 for (
const std::string &multname : mults) {
630 var_description &multdescr = variables.find(multname)->second;
631 dal::bit_vector kept;
632 size_type nbdof = multdescr.mf->nb_dof();
633 for (
const size_type &icol : glob_columns)
634 if (icol >= s && icol < s + nbdof) kept.add(icol-s);
635 if (multdescr.filter & VDESCRFILTER_REGION)
636 kept &= multdescr.mf->dof_on_region(multdescr.filter_region);
637 multdescr.partial_mf->adapt(kept);
638 multdescr.set_size();
639 multdescr.v_num = act_counter();
640 s += multdescr.mf->nb_dof();
648 resize_global_system();
660 if (variables.size() == 0)
661 ost <<
"Model with no variable nor data" << endl;
663 ost <<
"List of model variables and data:" << endl;
664 for (
int vartype=0; vartype < 3; ++vartype)
665 for (
const auto &v : variables) {
666 const var_description &vdescr = v.second;
667 bool is_variable = vdescr.is_variable;
670 if (!is_variable || is_disabled)
continue;
671 }
else if (vartype == 1) {
672 if (!is_disabled)
continue;
673 }
else if (vartype == 2) {
674 if (is_variable)
continue;
676 ost << (is_variable ?
"Variable " :
"Data ");
677 ost << std::setw(30) << std::left << v.first;
678 ost << std::setw(2) << std::right << vdescr.n_iter;
679 ost << ((vdescr.n_iter == 1) ?
" copy " :
" copies ");
680 ost << (vdescr.is_fem_dofs ?
"fem dependant " :
"constant size ");
681 ost << std::setw(8) << std::right << vdescr.size();
683 ost << ((vdescr.size() > 1) ?
" doubles." :
" double.");
684 ost << (is_disabled ?
"\t (disabled)" :
"\t ");
685 if (vdescr.imd != 0) ost <<
"\t (is im_data)";
686 if (vdescr.is_affine_dependent) ost <<
"\t (is affine dependent)";
689 for (
const auto &vargroup : variable_groups) {
690 ost <<
"Variable group " << std::setw(30) << std::left
692 if (vargroup.second.size()) {
694 for (
const std::string &vname : vargroup.second) {
695 ost << (first ?
" " :
", ") << vname;
700 ost <<
" empty" << endl;
705 void model::listresiduals(std::ostream &ost)
const {
707 if (variables.size() == 0)
708 ost <<
"Model with no variable nor data" << endl;
711 for (
const auto &v : variables) {
712 if (v.second.is_variable) {
713 const model_real_plain_vector &rhs = v.second.is_internal
715 const gmm::sub_interval &II = interval_of_variable(v.first);
716 scalar_type res = gmm::vect_norm2(gmm::sub_vector(rhs, II));
717 if (!firstvar) cout <<
", ";
718 ost <<
"res_" << v.first <<
"= " << std::setw(11) << res;
728 bgeot::multi_index sizes(1);
734 const bgeot::multi_index &sizes,
736 check_name_validity(name);
737 variables.emplace(name, var_description(
true,
is_complex(), 0, 0, niter));
738 variables[name].qdims = sizes;
739 act_size_to_be_done =
true;
740 variables[name].set_size();
741 GMM_ASSERT1(variables[name].qdim(),
742 "Variables of null size are not allowed");
747 bgeot::multi_index sizes(1);
753 const bgeot::multi_index &sizes) {
754 GMM_ASSERT1(!(variables[name].is_fem_dofs),
755 "Cannot explicitly resize a fem variable or data");
756 GMM_ASSERT1(variables[name].imd == 0,
757 "Cannot explicitly resize an im data");
758 variables[name].qdims = sizes;
759 variables[name].set_size();
764 bgeot::multi_index sizes(1);
770 const bgeot::multi_index &sizes,
772 check_name_validity(name);
773 variables.emplace(name, var_description(
false,
is_complex(), 0, 0, niter));
774 variables[name].qdims = sizes;
775 GMM_ASSERT1(variables[name].qdim(),
"Data of null size are not allowed");
776 variables[name].set_size();
780 const base_matrix &M) {
783 GMM_ASSERT1(!(
is_complex()),
"Sorry, complex version to be done");
788 const base_complex_matrix &M) {
791 GMM_ASSERT1(!(
is_complex()),
"Sorry, complex version to be done");
796 const base_tensor &t) {
798 GMM_ASSERT1(!(
is_complex()),
"Sorry, complex version to be done");
803 const base_complex_tensor &t) {
805 GMM_ASSERT1(!(
is_complex()),
"Sorry, complex version to be done");
811 check_name_validity(name);
812 variables.emplace(name,
813 var_description(
true,
is_complex(), 0, &imd, niter));
814 variables[name].set_size();
816 act_size_to_be_done =
true;
822 variables[name].is_internal =
true;
827 check_name_validity(name);
828 variables.emplace(name,
829 var_description(
false,
is_complex(), 0, &imd, niter));
830 variables[name].set_size();
836 check_name_validity(name);
837 variables.emplace(name, var_description(
true,
is_complex(), &mf, 0, niter,
839 variables[name].set_size();
841 act_size_to_be_done =
true;
842 leading_dim = std::max(leading_dim, mf.
linked_mesh().dim());
848 check_name_validity(name);
849 variables.emplace(name, var_description(
true,
is_complex(), &mf, 0, niter,
850 VDESCRFILTER_REGION, region));
851 variables[name].set_size();
852 act_size_to_be_done =
true;
857 const std::string &org_name,
859 check_name_validity(name);
860 VAR_SET::const_iterator it = find_variable(org_name);
861 GMM_ASSERT1(it->second.is_variable && !(it->second.is_affine_dependent),
862 "The original variable should be a variable");
863 variables.emplace(name, variables[org_name]);
864 variables[name].is_affine_dependent =
true;
865 variables[name].org_name = org_name;
866 variables[name].alpha = alpha;
867 variables[name].set_size();
872 bgeot::multi_index sizes(1);
878 const bgeot::multi_index &sizes,
size_type niter) {
879 check_name_validity(name);
880 variables.emplace(name, var_description(
false,
is_complex(), &mf, 0, niter,
882 variables[name].qdims = sizes;
883 GMM_ASSERT1(variables[name].qdim(),
"Data of null size are not allowed");
884 variables[name].set_size();
889 const std::string &primal_name,
891 check_name_validity(name);
892 variables.emplace(name, var_description(
true,
is_complex(), &mf, 0, niter,
895 variables[name].set_size();
896 act_size_to_be_done =
true;
901 size_type region,
const std::string &primal_name,
903 check_name_validity(name);
904 variables.emplace(name, var_description(
true,
is_complex(), &mf, 0, niter,
905 VDESCRFILTER_REGION_CTERM, region,
907 variables[name].set_size();
908 act_size_to_be_done =
true;
913 const std::string &primal_name,
916 check_name_validity(name);
917 variables.emplace(name, var_description(
true,
is_complex(), &mf, 0, niter,
918 VDESCRFILTER_INFSUP, region,
920 variables[name].set_size();
921 act_size_to_be_done =
true;
931 VAR_SET::iterator it = variables.find(name);
932 GMM_ASSERT1(it != variables.end(),
"Undefined variable " << name);
933 it->second.is_disabled = !enabled;
934 for (
auto &&v : variables) {
935 if (((v.second.filter & VDESCRFILTER_INFSUP) ||
936 (v.second.filter & VDESCRFILTER_CTERM))
937 && name.compare(v.second.filter_var) == 0) {
938 v.second.is_disabled = !enabled;
940 if (v.second.is_variable && v.second.is_affine_dependent
941 && name.compare(v.second.org_name) == 0)
942 v.second.is_disabled = !enabled;
944 if (!act_size_to_be_done) resize_global_system();
952 check_name_validity(name.substr(0, name.find(
"(")));
953 macro_dict.add_macro(name, expr);
957 { macro_dict.del_macro(name); }
960 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
961 valid_bricks.del(ib);
962 active_bricks.del(ib);
964 for (
size_type i = 0; i < bricks[ib].mims.size(); ++i) {
965 const mesh_im *mim = bricks[ib].mims[i];
967 for (dal::bv_visitor ibb(valid_bricks); !ibb.finished(); ++ibb) {
968 for (
size_type j = 0; j < bricks[ibb].mims.size(); ++j)
969 if (bricks[ibb].mims[j] == mim) found =
true;
971 for (
const auto &v : variables) {
972 if (v.second.is_fem_dofs &&
973 (v.second.filter & VDESCRFILTER_INFSUP) &&
974 mim == v.second.filter_mim) found =
true;
976 if (!found) sup_dependency(*mim);
979 is_linear_ = is_symmetric_ = is_coercive_ =
true;
980 for (dal::bv_visitor ibb(valid_bricks); !ibb.finished(); ++ibb) {
981 is_linear_ = is_linear_ && bricks[ibb].pbr->is_linear();
982 is_symmetric_ = is_symmetric_ && bricks[ibb].pbr->is_symmetric();
983 is_coercive_ = is_coercive_ && bricks[ibb].pbr->is_coercive();
985 bricks[ib] = brick_description();
989 for (dal::bv_visitor ibb(valid_bricks); !ibb.finished(); ++ibb) {
990 for (
const auto &vname : bricks[ibb].vlist)
991 GMM_ASSERT1(varname.compare(vname),
992 "Cannot delete a variable which is still used by a brick");
993 for (
const auto &dname : bricks[ibb].dlist)
994 GMM_ASSERT1(varname.compare(dname),
995 "Cannot delete a data which is still used by a brick");
998 VAR_SET::const_iterator it = find_variable(varname);
1000 if (it->second.is_fem_dofs) {
1001 const mesh_fem *mf = it->second.mf;
1003 for(VAR_SET::iterator it2 = variables.begin();
1004 it2 != variables.end(); ++it2) {
1005 if (it != it2 && it2->second.is_fem_dofs && mf == it2->second.mf)
1008 if (!found) sup_dependency(*mf);
1010 if (it->second.filter & VDESCRFILTER_INFSUP) {
1011 const mesh_im *mim = it->second.filter_mim;
1013 for (dal::bv_visitor ibb(valid_bricks); !ibb.finished(); ++ibb) {
1014 for (
size_type j = 0; j < bricks[ibb].mims.size(); ++j)
1015 if (bricks[ibb].mims[j] == mim) found =
true;
1017 for (VAR_SET::iterator it2 = variables.begin();
1018 it2 != variables.end(); ++it2) {
1019 if (it != it2 && it2->second.is_fem_dofs &&
1020 (it2->second.filter & VDESCRFILTER_INFSUP) &&
1021 mim == it2->second.filter_mim) found =
true;
1023 if (!found) sup_dependency(*mim);
1027 if (it->second.imd != 0) sup_dependency(*(it->second.imd));
1029 variables.erase(varname);
1030 act_size_to_be_done =
true;
1034 const varnamelist &datanames,
1035 const termlist &terms,
1036 const mimlist &mims,
size_type region) {
1037 size_type ib = valid_bricks.first_false();
1039 for (
size_type i = 0; i < terms.size(); ++i)
1040 if (terms[i].is_global && terms[i].is_matrix_term && pbr->is_linear())
1041 GMM_ASSERT1(
false,
"Global linear matrix terms are not allowed");
1043 if (ib == bricks.size())
1044 bricks.push_back(brick_description(pbr, varnames, datanames, terms,
1047 bricks[ib] = brick_description(pbr, varnames, datanames, terms,
1049 active_bricks.add(ib);
1050 valid_bricks.add(ib);
1057 "Impossible to add a complex brick to a real model");
1059 bricks[ib].cmatlist.resize(terms.size());
1060 bricks[ib].cveclist[0].resize(terms.size());
1061 bricks[ib].cveclist_sym[0].resize(terms.size());
1063 bricks[ib].rmatlist.resize(terms.size());
1064 bricks[ib].rveclist[0].resize(terms.size());
1065 bricks[ib].rveclist_sym[0].resize(terms.size());
1067 is_linear_ = is_linear_ && pbr->is_linear();
1068 is_symmetric_ = is_symmetric_ && pbr->is_symmetric();
1069 is_coercive_ = is_coercive_ && pbr->is_coercive();
1071 for (
const auto &vname : varnames)
1072 GMM_ASSERT1(variables.count(vname),
1073 "Undefined model variable " << vname);
1075 for (
const auto &dname : datanames)
1076 GMM_ASSERT1(variables.count(dname),
1077 "Undefined model data or variable " << dname);
1083 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
1085 bricks[ib].mims.push_back(&mim);
1086 add_dependency(mim);
1090 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
1092 bricks[ib].tlist = terms;
1093 if (
is_complex() && bricks[ib].pbr->is_complex()) {
1094 bricks.back().cmatlist.resize(terms.size());
1095 bricks.back().cveclist[0].resize(terms.size());
1096 bricks.back().cveclist_sym[0].resize(terms.size());
1098 bricks.back().rmatlist.resize(terms.size());
1099 bricks.back().rveclist[0].resize(terms.size());
1100 bricks.back().rveclist_sym[0].resize(terms.size());
1105 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
1107 bricks[ib].vlist = vl;
1108 for (
const auto &v : vl)
1109 GMM_ASSERT1(variables.count(v),
"Undefined model variable " << v);
1113 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
1115 bricks[ib].dlist = dl;
1116 for (
const auto &v : dl)
1117 GMM_ASSERT1(variables.count(v),
"Undefined model variable " << v);
1121 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
1123 bricks[ib].mims = ml;
1124 for (
const auto &mim : ml) add_dependency(*mim);
1128 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
1130 bricks[ib].is_update_brick = flag;
1133 void model::set_time(scalar_type t,
bool to_init) {
1134 static const std::string varname(
"t");
1135 VAR_SET::iterator it = variables.find(varname);
1136 if (it == variables.end()) {
1139 GMM_ASSERT1(it->second.size() == 1,
"Time data should be of size 1");
1141 if (it == variables.end() || to_init) {
1149 scalar_type model::get_time() {
1150 static const std::string varname(
"t");
1151 set_time(scalar_type(0),
false);
1158 void model::call_init_affine_dependent_variables(
int version) {
1159 for (VAR_SET::iterator it = variables.begin();
1160 it != variables.end(); ++it) {
1161 var_description &vdescr = it->second;
1162 if (vdescr.is_variable && vdescr.ptsc) {
1164 vdescr.ptsc->init_affine_dependent_variables_precomputation(*
this);
1166 vdescr.ptsc->init_affine_dependent_variables(*
this);
1171 void model::shift_variables_for_time_integration() {
1172 for (VAR_SET::iterator it = variables.begin();
1173 it != variables.end(); ++it)
1174 if (it->second.is_variable && it->second.ptsc)
1175 it->second.ptsc->shift_variables(*
this);
1178 void model::add_time_integration_scheme(
const std::string &varname,
1179 ptime_scheme ptsc) {
1180 VAR_SET::iterator it = variables.find(varname);
1181 GMM_ASSERT1(it != variables.end(),
"Undefined variable " << varname);
1182 GMM_ASSERT1(it->second.is_variable && !(it->second.is_affine_dependent),
1183 "Cannot apply an integration scheme to " << varname);
1184 it->second.ptsc = ptsc;
1191 time_integration = 1;
1194 void model::copy_init_time_derivative() {
1196 for (VAR_SET::iterator it = variables.begin();
1197 it != variables.end(); ++it)
1198 if (it->second.is_variable && it->second.ptsc) {
1200 std::string name_v, name_previous_v;
1201 it->second.ptsc->time_derivative_to_be_initialized(name_v,
1204 if (name_v.size()) {
1222 class APIDECL first_order_theta_method_scheme
1223 :
public virtual_time_scheme {
1225 std::string U, U0, V, V0;
1230 virtual void init_affine_dependent_variables(model &md)
const {
1231 scalar_type dt = md.get_time_step();
1232 scalar_type a = scalar_type(1)/(theta*dt);
1233 scalar_type b = (scalar_type(1)-theta)/theta;
1234 md.set_factor_of_variable(V, a);
1235 if (md.is_complex()) {
1236 gmm::add(gmm::scaled(md.complex_variable(U0), -complex_type(a)),
1237 gmm::scaled(md.complex_variable(V0), -complex_type(b)),
1238 md.set_complex_constant_part(V));
1241 gmm::add(gmm::scaled(md.real_variable(U0), -a),
1242 gmm::scaled(md.real_variable(V0), -b),
1243 md.set_real_constant_part(V));
1248 virtual void init_affine_dependent_variables_precomputation(model &md)
1250 scalar_type dt = md.get_time_step();
1251 md.set_factor_of_variable(V, scalar_type(1)/dt);
1252 if (md.is_complex()) {
1253 gmm::copy(gmm::scaled(md.complex_variable(U0), -complex_type(1)/dt),
1254 md.set_complex_constant_part(V));
1257 gmm::copy(gmm::scaled(md.real_variable(U0), -scalar_type(1)/dt),
1258 md.set_real_constant_part(V));
1262 virtual void time_derivative_to_be_initialized
1263 (std::string &name_v, std::string &name_previous_v)
const
1264 {
if (theta != scalar_type(1)) { name_v = V; name_previous_v = V0; } }
1266 virtual void shift_variables(model &md)
const {
1267 if (md.is_complex()) {
1268 gmm::copy(md.complex_variable(U), md.set_complex_variable(U0));
1269 gmm::copy(md.complex_variable(V), md.set_complex_variable(V0));
1271 gmm::copy(md.real_variable(U), md.set_real_variable(U0));
1272 gmm::copy(md.real_variable(V), md.set_real_variable(V0));
1277 first_order_theta_method_scheme(model &md, std::string varname,
1280 U0 =
"Previous_" + U;
1282 V0 =
"Previous_Dot_" + U;
1284 GMM_ASSERT1(theta > scalar_type(0) && theta <= scalar_type(1),
1285 "Invalid value of theta parameter for the theta-method");
1287 if (!(md.variable_exists(V)))
1288 md.add_affine_dependent_variable(V, U);
1289 const mesh_fem *mf = md.pmesh_fem_of_variable(U);
1290 size_type s = md.is_complex() ? gmm::vect_size(md.complex_variable(U))
1291 : gmm::vect_size(md.real_variable(U));
1294 if (!(md.variable_exists(U0))) md.add_fem_data(U0, *mf);
1295 if (!(md.variable_exists(V0))) md.add_fem_data(V0, *mf);
1297 if (!(md.variable_exists(U0))) md.add_fixed_size_data(U0, s);
1298 if (!(md.variable_exists(V0))) md.add_fixed_size_data(V0, s);
1305 void add_theta_method_for_first_order(model &md,
const std::string &varname,
1306 scalar_type theta) {
1308 = std::make_shared<first_order_theta_method_scheme>(md, varname,theta);
1309 md.add_time_integration_scheme(varname, ptsc);
1318 class APIDECL second_order_theta_method_scheme
1319 :
public virtual_time_scheme {
1321 std::string U, U0, V, V0, A, A0;
1327 virtual void init_affine_dependent_variables(model &md)
const {
1328 scalar_type dt = md.get_time_step();
1329 md.set_factor_of_variable(V, scalar_type(1)/(theta*dt));
1330 md.set_factor_of_variable(A, scalar_type(1)/(theta*theta*dt*dt));
1331 if (md.is_complex()) {
1332 gmm::add(gmm::scaled(md.complex_variable(U0),
1333 -complex_type(1)/(theta*dt)),
1334 gmm::scaled(md.complex_variable(V0),
1335 -(complex_type(1)-complex_type(theta))/theta),
1336 md.set_complex_constant_part(V));
1337 gmm::add(gmm::scaled(md.complex_variable(U0),
1338 -complex_type(1)/(theta*theta*dt*dt)),
1339 gmm::scaled(md.complex_variable(A0),
1340 -(complex_type(1)-complex_type(theta))/theta),
1341 md.set_complex_constant_part(A));
1342 gmm::add(gmm::scaled(md.complex_variable(V0),
1343 -complex_type(1)/(theta*theta*dt)),
1344 md.set_complex_constant_part(A));
1348 gmm::add(gmm::scaled(md.real_variable(U0),
1349 -scalar_type(1)/(theta*dt)),
1350 gmm::scaled(md.real_variable(V0),
1351 -(scalar_type(1)-theta)/theta),
1352 md.set_real_constant_part(V));
1353 gmm::add(gmm::scaled(md.real_variable(U0),
1354 -scalar_type(1)/(theta*theta*dt*dt)),
1355 gmm::scaled(md.real_variable(A0),
1356 -(scalar_type(1)-theta)/theta),
1357 md.set_real_constant_part(A));
1358 gmm::add(gmm::scaled(md.real_variable(V0),
1359 -scalar_type(1)/(theta*theta*dt)),
1360 md.set_real_constant_part(A));
1367 virtual void init_affine_dependent_variables_precomputation(model &md)
1369 scalar_type dt = md.get_time_step();
1370 md.set_factor_of_variable(V, scalar_type(1)/dt);
1371 md.set_factor_of_variable(A, scalar_type(1)/(dt*dt));
1372 if (md.is_complex()) {
1373 gmm::copy(gmm::scaled(md.complex_variable(U0),
1374 -complex_type(1)/dt),
1375 md.set_complex_constant_part(V));
1376 gmm::add(gmm::scaled(md.complex_variable(U0),
1377 -complex_type(1)/(dt*dt)),
1378 gmm::scaled(md.complex_variable(V0),
1379 -complex_type(1)/dt),
1380 md.set_complex_constant_part(A));
1382 gmm::copy(gmm::scaled(md.real_variable(U0),
1383 -scalar_type(1)/dt),
1384 md.set_real_constant_part(V));
1385 gmm::add(gmm::scaled(md.real_variable(U0),
1386 -scalar_type(1)/(dt*dt)),
1387 gmm::scaled(md.real_variable(V0),
1388 -scalar_type(1)/dt),
1389 md.set_real_constant_part(A));
1393 virtual void time_derivative_to_be_initialized
1394 (std::string &name_v, std::string &name_previous_v)
const
1395 {
if (theta != scalar_type(1)) { name_v = A; name_previous_v = A0; } }
1397 virtual void shift_variables(model &md)
const {
1398 if (md.is_complex()) {
1399 gmm::copy(md.complex_variable(U), md.set_complex_variable(U0));
1400 gmm::copy(md.complex_variable(V), md.set_complex_variable(V0));
1401 gmm::copy(md.complex_variable(A), md.set_complex_variable(A0));
1403 gmm::copy(md.real_variable(U), md.set_real_variable(U0));
1404 gmm::copy(md.real_variable(V), md.set_real_variable(V0));
1405 gmm::copy(md.real_variable(A), md.set_real_variable(A0));
1410 second_order_theta_method_scheme(model &md, std::string varname,
1413 U0 =
"Previous_" + U;
1415 V0 =
"Previous_Dot_" + U;
1417 A0 =
"Previous_Dot2_" + U;
1419 GMM_ASSERT1(theta > scalar_type(0) && theta <= scalar_type(1),
1420 "Invalid value of theta parameter for the theta-method");
1422 if (!(md.variable_exists(V)))
1423 md.add_affine_dependent_variable(V, U);
1424 if (!(md.variable_exists(A)))
1425 md.add_affine_dependent_variable(A, U);
1427 const mesh_fem *mf = md.pmesh_fem_of_variable(U);
1428 size_type s = md.is_complex() ? gmm::vect_size(md.complex_variable(U))
1429 : gmm::vect_size(md.real_variable(U));
1432 if (!(md.variable_exists(U0))) md.add_fem_data(U0, *mf);
1433 if (!(md.variable_exists(V0))) md.add_fem_data(V0, *mf);
1434 if (!(md.variable_exists(A0))) md.add_fem_data(A0, *mf);
1436 if (!(md.variable_exists(U0))) md.add_fixed_size_data(U0, s);
1437 if (!(md.variable_exists(V0))) md.add_fixed_size_data(V0, s);
1438 if (!(md.variable_exists(A0))) md.add_fixed_size_data(A0, s);
1445 void add_theta_method_for_second_order(model &md,
const std::string &varname,
1446 scalar_type theta) {
1447 ptime_scheme ptsc = std::make_shared<second_order_theta_method_scheme>
1449 md.add_time_integration_scheme(varname, ptsc);
1459 class APIDECL Newmark_scheme
1460 :
public virtual_time_scheme {
1462 std::string U, U0, V, V0, A, A0;
1463 scalar_type beta, gamma;
1468 virtual void init_affine_dependent_variables(model &md)
const {
1469 scalar_type dt = md.get_time_step();
1470 scalar_type a0 = scalar_type(1)/(beta*dt*dt), a1 = dt*a0;
1471 scalar_type a2 = (scalar_type(1) - scalar_type(2)*beta)
1472 / (scalar_type(2)*beta);
1473 scalar_type b0 = gamma/(beta*dt), b1 = (beta-gamma)/beta;
1474 scalar_type b2 = dt*(scalar_type(1)-gamma/(scalar_type(2)*beta));
1476 md.set_factor_of_variable(V, b0);
1477 md.set_factor_of_variable(A, a0);
1478 if (md.is_complex()) {
1479 gmm::add(gmm::scaled(md.complex_variable(U0), -complex_type(b0)),
1480 gmm::scaled(md.complex_variable(V0), complex_type(b1)),
1481 md.set_complex_constant_part(V));
1482 gmm::add(gmm::scaled(md.complex_variable(A0), complex_type(b2)),
1483 md.set_complex_constant_part(V));
1484 gmm::add(gmm::scaled(md.complex_variable(U0), -complex_type(a0)),
1485 gmm::scaled(md.complex_variable(V0), -complex_type(a1)),
1486 md.set_complex_constant_part(A));
1487 gmm::add(gmm::scaled(md.complex_variable(A0), -complex_type(a2)),
1488 md.set_complex_constant_part(A));
1490 gmm::add(gmm::scaled(md.real_variable(U0), -b0),
1491 gmm::scaled(md.real_variable(V0), b1),
1492 md.set_real_constant_part(V));
1493 gmm::add(gmm::scaled(md.real_variable(A0), b2),
1494 md.set_real_constant_part(V));
1495 gmm::add(gmm::scaled(md.real_variable(U0), -a0),
1496 gmm::scaled(md.real_variable(V0), -a1),
1497 md.set_real_constant_part(A));
1498 gmm::add(gmm::scaled(md.real_variable(A0), -a2),
1499 md.set_real_constant_part(A));
1506 virtual void init_affine_dependent_variables_precomputation(model &md)
1508 scalar_type dt = md.get_time_step();
1509 md.set_factor_of_variable(V, scalar_type(1)/dt);
1510 md.set_factor_of_variable(A, scalar_type(1)/(dt*dt));
1511 if (md.is_complex()) {
1512 gmm::copy(gmm::scaled(md.complex_variable(U0),
1513 -complex_type(1)/dt),
1514 md.set_complex_constant_part(V));
1515 gmm::add(gmm::scaled(md.complex_variable(U0),
1516 -complex_type(1)/(dt*dt)),
1517 gmm::scaled(md.complex_variable(V0),
1518 -complex_type(1)/dt),
1519 md.set_complex_constant_part(A));
1521 gmm::copy(gmm::scaled(md.real_variable(U0),
1522 -scalar_type(1)/dt),
1523 md.set_real_constant_part(V));
1524 gmm::add(gmm::scaled(md.real_variable(U0),
1525 -scalar_type(1)/(dt*dt)),
1526 gmm::scaled(md.real_variable(V0),
1527 -scalar_type(1)/dt),
1528 md.set_real_constant_part(A));
1532 virtual void time_derivative_to_be_initialized
1533 (std::string &name_v, std::string &name_previous_v)
const {
1534 if (beta != scalar_type(0.5) || gamma != scalar_type(1))
1535 { name_v = A; name_previous_v = A0; }
1538 virtual void shift_variables(model &md)
const {
1539 if (md.is_complex()) {
1540 gmm::copy(md.complex_variable(U), md.set_complex_variable(U0));
1541 gmm::copy(md.complex_variable(V), md.set_complex_variable(V0));
1542 gmm::copy(md.complex_variable(A), md.set_complex_variable(A0));
1544 gmm::copy(md.real_variable(U), md.set_real_variable(U0));
1545 gmm::copy(md.real_variable(V), md.set_real_variable(V0));
1546 gmm::copy(md.real_variable(A), md.set_real_variable(A0));
1551 Newmark_scheme(model &md, std::string varname,
1552 scalar_type be, scalar_type ga) {
1554 U0 =
"Previous_" + U;
1556 V0 =
"Previous_Dot_" + U;
1558 A0 =
"Previous_Dot2_" + U;
1559 beta = be; gamma = ga;
1560 GMM_ASSERT1(beta > scalar_type(0) && beta <= scalar_type(1)
1561 && gamma >= scalar_type(0.5) && gamma <= scalar_type(1),
1562 "Invalid parameter values for the Newmark scheme");
1564 if (!(md.variable_exists(V)))
1565 md.add_affine_dependent_variable(V, U);
1566 if (!(md.variable_exists(A)))
1567 md.add_affine_dependent_variable(A, U);
1569 const mesh_fem *mf = md.pmesh_fem_of_variable(U);
1570 size_type s = md.is_complex() ? gmm::vect_size(md.complex_variable(U))
1571 : gmm::vect_size(md.real_variable(U));
1574 if (!(md.variable_exists(U0))) md.add_fem_data(U0, *mf);
1575 if (!(md.variable_exists(V0))) md.add_fem_data(V0, *mf);
1576 if (!(md.variable_exists(A0))) md.add_fem_data(A0, *mf);
1578 if (!(md.variable_exists(U0))) md.add_fixed_size_data(U0, s);
1579 if (!(md.variable_exists(V0))) md.add_fixed_size_data(V0, s);
1580 if (!(md.variable_exists(A0))) md.add_fixed_size_data(A0, s);
1587 void add_Newmark_scheme(model &md,
const std::string &varname,
1588 scalar_type beta, scalar_type gamma) {
1589 ptime_scheme ptsc = std::make_shared<Newmark_scheme>
1590 (md, varname, beta, gamma);
1591 md.add_time_integration_scheme(varname, ptsc);
1600 class APIDECL Houbolt_scheme
1601 :
public virtual_time_scheme {
1603 std::string U, U01, U02, U03, V, A;
1608 virtual void init_affine_dependent_variables(model &md)
const {
1609 scalar_type dt = md.get_time_step();
1610 scalar_type a0 = scalar_type(2)/(dt*dt);
1611 scalar_type a1 = scalar_type(5)/(dt*dt);
1612 scalar_type a2 = scalar_type(4)/(dt*dt);
1613 scalar_type a3 = scalar_type(1)/(dt*dt);
1614 scalar_type b0 = scalar_type(11)/(scalar_type(6)*dt);
1615 scalar_type b1 = scalar_type(18)/(scalar_type(6)*dt);
1616 scalar_type b2 = scalar_type(9)/(scalar_type(6)*dt);
1617 scalar_type b3 = scalar_type(2)/(scalar_type(6)*dt);
1619 md.set_factor_of_variable(V, b0);
1620 md.set_factor_of_variable(A, a0);
1621 if (md.is_complex()) {
1622 gmm::add(gmm::scaled(md.complex_variable(U01), -complex_type(b1)),
1623 gmm::scaled(md.complex_variable(U02), complex_type(b2)),
1624 md.set_complex_constant_part(V));
1625 gmm::add(gmm::scaled(md.complex_variable(U03), -complex_type(b3)),
1626 md.set_complex_constant_part(V));
1627 gmm::add(gmm::scaled(md.complex_variable(U01), -complex_type(a1)),
1628 gmm::scaled(md.complex_variable(U02), complex_type(a2)),
1629 md.set_complex_constant_part(A));
1630 gmm::add(gmm::scaled(md.complex_variable(U03), -complex_type(a3)),
1631 md.set_complex_constant_part(A));
1633 gmm::add(gmm::scaled(md.real_variable(U01), -b1),
1634 gmm::scaled(md.real_variable(U02), b2),
1635 md.set_real_constant_part(V));
1636 gmm::add(gmm::scaled(md.real_variable(U03), -b3),
1637 md.set_real_constant_part(V));
1638 gmm::add(gmm::scaled(md.real_variable(U01), -a1),
1639 gmm::scaled(md.real_variable(U02), a2),
1640 md.set_real_constant_part(A));
1641 gmm::add(gmm::scaled(md.real_variable(U03), -a3),
1642 md.set_real_constant_part(A));
1646 virtual void init_affine_dependent_variables_precomputation(model &md)
1651 virtual void time_derivative_to_be_initialized
1652 (std::string &name_v, std::string &name_previous_v)
const {
1654 (void) name_previous_v;
1657 virtual void shift_variables(model &md)
const {
1658 if (md.is_complex()) {
1659 gmm::copy(md.complex_variable(U02), md.set_complex_variable(U03));
1660 gmm::copy(md.complex_variable(U01), md.set_complex_variable(U02));
1661 gmm::copy(md.complex_variable(U), md.set_complex_variable(U01));
1663 gmm::copy(md.real_variable(U02), md.set_real_variable(U03));
1664 gmm::copy(md.real_variable(U01), md.set_real_variable(U02));
1665 gmm::copy(md.real_variable(U), md.set_real_variable(U01));
1670 Houbolt_scheme(model &md, std::string varname) {
1672 U01 =
"Previous_" + U;
1673 U02 =
"Previous2_" + U;
1674 U03 =
"Previous3_" + U;
1678 if (!(md.variable_exists(V)))
1679 md.add_affine_dependent_variable(V, U);
1680 if (!(md.variable_exists(A)))
1681 md.add_affine_dependent_variable(A, U);
1683 const mesh_fem *mf = md.pmesh_fem_of_variable(U);
1684 size_type s = md.is_complex() ? gmm::vect_size(md.complex_variable(U))
1685 : gmm::vect_size(md.real_variable(U));
1688 if (!(md.variable_exists(U01))) md.add_fem_data(U01, *mf);
1689 if (!(md.variable_exists(U02))) md.add_fem_data(U02, *mf);
1690 if (!(md.variable_exists(U03))) md.add_fem_data(U03, *mf);
1692 if (!(md.variable_exists(U01))) md.add_fixed_size_data(U01, s);
1693 if (!(md.variable_exists(U02))) md.add_fixed_size_data(U02, s);
1694 if (!(md.variable_exists(U03))) md.add_fixed_size_data(U03, s);
1701 void add_Houbolt_scheme(model &md,
const std::string &varname) {
1702 ptime_scheme ptsc = std::make_shared<Houbolt_scheme>
1704 md.add_time_integration_scheme(varname, ptsc);
1716 GMM_ASSERT1(valid_bricks[ibrick],
"Inexistent brick");
1717 pbrick pbr = bricks[ibrick].pbr;
1719 bricks[ibrick].pdispatch = pdispatch;
1722 = std::max(
size_type(1), pdispatch->nbrhs());
1724 gmm::resize(bricks[ibrick].coeffs, nbrhs);
1727 bricks[ibrick].cveclist.resize(nbrhs);
1728 bricks[ibrick].cveclist_sym.resize(nbrhs);
1730 bricks[ibrick].cveclist[k] = bricks[ibrick].cveclist[0];
1731 bricks[ibrick].cveclist_sym[k] = bricks[ibrick].cveclist_sym[0];
1734 bricks[ibrick].rveclist.resize(nbrhs);
1735 bricks[ibrick].rveclist_sym.resize(nbrhs);
1737 bricks[ibrick].rveclist[k] = bricks[ibrick].rveclist[0];
1738 bricks[ibrick].rveclist_sym[k] = bricks[ibrick].rveclist_sym[0];
1745 GMM_ASSERT1(valid_bricks[ind_brick],
"Inexistent brick");
1746 GMM_ASSERT1(ind_var < bricks[ind_brick].vlist.size(),
1747 "Inexistent brick variable");
1748 return bricks[ind_brick].vlist[ind_var];
1753 GMM_ASSERT1(valid_bricks[ind_brick],
"Inexistent brick");
1754 GMM_ASSERT1(ind_data < bricks[ind_brick].dlist.size(),
1755 "Inexistent brick data");
1756 return bricks[ind_brick].dlist[ind_data];
1760 if (valid_bricks.card() == 0)
1761 ost <<
"Model with no bricks" << endl;
1763 ost <<
"List of model bricks:" << endl;
1764 for (dal::bv_visitor i(valid_bricks); !i.finished(); ++i) {
1765 ost <<
"Brick " << std::setw(3) << std::right << i + base_id
1766 <<
" " << std::setw(20) << std::right
1767 << bricks[i].pbr->brick_name();
1768 if (!(active_bricks[i])) ost <<
" (deactivated)";
1769 if (bricks[i].pdispatch) ost <<
" (dispatched)";
1770 ost << endl <<
" concerned variables: " << bricks[i].vlist[0];
1771 for (
size_type j = 1; j < bricks[i].vlist.size(); ++j)
1772 ost <<
", " << bricks[i].vlist[j];
1774 ost <<
" brick with " << bricks[i].tlist.size() <<
" term";
1775 if (bricks[i].tlist.size() > 1) ost <<
"s";
1784 void model::brick_init(
size_type ib, build_version version,
1786 const brick_description &brick = bricks[ib];
1787 bool cplx =
is_complex() && brick.pbr->is_complex();
1790 for (
size_type j = 0; j < brick.tlist.size(); ++j) {
1791 const term_description &term = brick.tlist[j];
1792 bool isg = term.is_global;
1794 gmm::vect_size(crhs) : gmm::vect_size(rrhs);
1795 size_type nbd1 = isg ? nbgdof : variables[term.var1].size();
1796 size_type nbd2 = isg ? nbgdof : (term.is_matrix_term ?
1797 variables[term.var2].size() : 0);
1798 if (term.is_matrix_term &&
1799 (brick.pbr->is_linear() || (version | BUILD_MATRIX))) {
1800 if (version | BUILD_ON_DATA_CHANGE) {
1802 gmm::resize(brick.cmatlist[j], nbd1, nbd2);
1804 gmm::resize(brick.rmatlist[j], nbd1, nbd2);
1807 brick.cmatlist[j] = model_complex_sparse_matrix(nbd1, nbd2);
1809 brick.rmatlist[j] = model_real_sparse_matrix(nbd1, nbd2);
1812 if (brick.pbr->is_linear() || (version | BUILD_RHS)) {
1813 for (
size_type k = 0; k < brick.nbrhs; ++k) {
1815 if (k == rhs_ind)
gmm::clear(brick.cveclist[k][j]);
1817 if (term.is_symmetric && term.var1.compare(term.var2)) {
1818 if (k == rhs_ind)
gmm::clear(brick.cveclist_sym[k][j]);
1822 if (k == rhs_ind)
gmm::clear(brick.rveclist[k][j]);
1824 if (term.is_symmetric && term.var1.compare(term.var2)) {
1825 if (k == rhs_ind)
gmm::clear(brick.rveclist_sym[k][j]);
1834 void model::post_to_variables_step(){}
1836 void model::brick_call(
size_type ib, build_version version,
1839 const brick_description &brick = bricks[ib];
1840 bool cplx =
is_complex() && brick.pbr->is_complex();
1842 brick_init(ib, version, rhs_ind);
1846 brick.pbr->complex_pre_assembly_in_serial(*
this, ib, brick.vlist,
1847 brick.dlist, brick.mims,
1849 brick.cveclist[rhs_ind],
1850 brick.cveclist_sym[rhs_ind],
1851 brick.region, version);
1856 accumulated_distro<complex_matlist> cmatlist(brick.cmatlist);
1857 accumulated_distro<complex_veclist> cveclist(brick.cveclist[rhs_ind]);
1858 accumulated_distro<complex_veclist> cveclist_sym(brick.cveclist_sym[rhs_ind]);
1862 brick.pbr->asm_complex_tangent_terms(*
this, ib, brick.vlist,
1863 brick.dlist, brick.mims,
1867 brick.region, version);
1870 brick.pbr->complex_post_assembly_in_serial(*
this, ib, brick.vlist,
1871 brick.dlist, brick.mims,
1873 brick.cveclist[rhs_ind],
1874 brick.cveclist_sym[rhs_ind],
1875 brick.region, version);
1877 if (brick.is_update_brick)
1879 for (
auto &&mat : brick.cmatlist)
1882 for (
auto &&vecs : brick.cveclist)
1883 for (
auto &&vec : vecs)
1886 for (
auto &&vecs : brick.cveclist_sym)
1887 for (
auto &&vec : vecs)
1893 brick.pbr->real_pre_assembly_in_serial(*
this, ib, brick.vlist,
1894 brick.dlist, brick.mims,
1896 brick.rveclist[rhs_ind],
1897 brick.rveclist_sym[rhs_ind],
1898 brick.region, version);
1901 accumulated_distro<real_matlist> rmatlist(brick.rmatlist);
1902 accumulated_distro<real_veclist> rveclist(brick.rveclist[rhs_ind]);
1903 accumulated_distro<real_veclist> rveclist_sym(brick.rveclist_sym[rhs_ind]);
1907 brick.pbr->asm_real_tangent_terms(*
this, ib, brick.vlist,
1908 brick.dlist, brick.mims,
1916 brick.pbr->real_post_assembly_in_serial(*
this, ib, brick.vlist,
1917 brick.dlist, brick.mims,
1919 brick.rveclist[rhs_ind],
1920 brick.rveclist_sym[rhs_ind],
1921 brick.region, version);
1923 if (brick.is_update_brick)
1925 for (
auto &&mat : brick.rmatlist)
1928 for (
auto &&vecs : brick.rveclist)
1929 for (
auto &&vec : vecs)
1932 for (
auto &&vecs : brick.rveclist_sym)
1933 for (
auto &&vec : vecs)
1940 void model::set_dispatch_coeff() {
1941 for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) {
1942 brick_description &brick = bricks[ib];
1943 if (brick.pdispatch)
1944 brick.pdispatch->set_dispatch_coeff(*
this, ib);
1950 context_check();
if (act_size_to_be_done) actualize_sizes();
1951 for (
auto && v : variables) v.second.clear_temporaries();
1953 set_dispatch_coeff();
1955 for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) {
1956 brick_description &brick = bricks[ib];
1957 if (brick.pdispatch) {
1959 brick.pdispatch->next_complex_iter(*
this, ib, brick.vlist,
1961 brick.cmatlist, brick.cveclist,
1962 brick.cveclist_sym,
true);
1964 brick.pdispatch->next_real_iter(*
this, ib, brick.vlist, brick.dlist,
1965 brick.rmatlist, brick.rveclist,
1966 brick.rveclist_sym,
true);
1972 context_check();
if (act_size_to_be_done) actualize_sizes();
1973 set_dispatch_coeff();
1975 for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) {
1976 brick_description &brick = bricks[ib];
1977 if (brick.pdispatch) {
1979 brick.pdispatch->next_complex_iter(*
this, ib, brick.vlist,
1981 brick.cmatlist, brick.cveclist,
1982 brick.cveclist_sym,
false);
1984 brick.pdispatch->next_real_iter(*
this, ib, brick.vlist, brick.dlist,
1985 brick.rmatlist, brick.rveclist,
1986 brick.rveclist_sym,
false);
1990 for (
auto &&v : variables)
1991 for (
size_type i = 1; i < v.second.n_iter; ++i) {
1993 gmm::copy(v.second.complex_value[i-1], v.second.complex_value[i]);
1995 gmm::copy(v.second.real_value[i-1], v.second.real_value[i]);
1996 v.second.v_num_data[i] = act_counter();
2000 bool model::is_var_newer_than_brick(
const std::string &varname,
2002 const brick_description &brick = bricks[ib];
2003 var_description &vd = variables[varname];
2004 if (niter ==
size_type(-1)) niter = vd.default_iter;
2005 return (vd.v_num > brick.v_num || vd.v_num_data[niter] > brick.v_num);
2008 bool model::is_var_mf_newer_than_brick(
const std::string &varname,
2010 const brick_description &brick = bricks[ib];
2011 var_description &vd = variables[varname];
2012 return (vd.v_num > brick.v_num);
2015 bool model::is_mim_newer_than_brick(
const mesh_im &im,
2017 const brick_description &brick = bricks[ib];
2018 return (im.version_number() > brick.v_num);
2021 void model::define_variable_group(
const std::string &group_name,
2022 const std::vector<std::string> &nl) {
2024 "variables cannot be the same as a variable name");
2026 std::set<const mesh *> ms;
2027 bool is_data_ =
false;
2028 for (
size_type i = 0; i < nl.size(); ++i) {
2033 "It is not possible to mix variables and data in a group");
2036 "All variables in a group have to exist in the model");
2038 GMM_ASSERT1(mf,
"Variables in a group should be fem variables");
2039 GMM_ASSERT1(ms.find(&(mf->linked_mesh())) == ms.end(),
2040 "Two variables in a group cannot share the same mesh");
2041 ms.insert(&(mf->linked_mesh()));
2043 variable_groups[group_name] = nl;
2046 void model::add_assembly_assignments(
const std::string &varname,
2049 GMM_ASSERT1(order < 3 || order ==
size_type(-1),
"Bad order value");
2050 const im_data *imd = pim_data_of_variable(varname);
2051 GMM_ASSERT1(imd != 0,
"Only applicable to im_data");
2052 assignement_desc as;
2053 as.varname = varname; as.expr = expr; as.region = rg; as.order = order;
2055 assignments.push_back(as);
2058 void model::add_temporaries(
const varnamelist &vl,
2059 gmm::uint64_type id_num)
const {
2060 for (
size_type i = 0; i < vl.size(); ++i) {
2061 var_description &vd = variables[vl[i]];
2062 if (vd.n_iter > 1) {
2063 vd.add_temporary(id_num);
2068 bool model::temporary_uptodate(
const std::string &varname,
2069 gmm::uint64_type id_num,
2071 var_description &vd = variables[varname];
2073 for (; ind < vd.n_iter + vd.n_temp_iter ; ++ind) {
2074 if (vd.v_num_var_iter[ind] == id_num)
break;
2076 if (ind < vd.n_iter + vd.n_temp_iter) {
2077 if (vd.v_num_iter[ind] <= vd.v_num_data[vd.default_iter]) {
2078 vd.v_num_iter[ind] = act_counter();
2087 void model::set_default_iter_of_variable(
const std::string &varname,
2090 var_description &vd = variables[varname];
2091 GMM_ASSERT1(ind < vd.n_iter + vd.n_temp_iter,
2092 "Inexistent iteration " << ind);
2093 vd.default_iter = ind;
2097 void model::reset_default_iter_of_variables(
const varnamelist &vl)
const {
2098 for (
size_type i = 0; i < vl.size(); ++i)
2099 variables[vl[i]].default_iter = 0;
2102 const model_real_sparse_matrix &
2104 GMM_ASSERT1(bricks[ib].tlist[iterm].is_matrix_term,
2105 "Not a matrix term !");
2106 GMM_ASSERT1(bricks[ib].pbr->is_linear(),
"Nonlinear term !");
2107 return bricks[ib].rmatlist[iterm];
2110 const model_complex_sparse_matrix &
2112 GMM_ASSERT1(bricks[ib].tlist[iterm].is_matrix_term,
2113 "Not a matrix term !");
2114 GMM_ASSERT1(bricks[ib].pbr->is_linear(),
"Nonlinear term !");
2115 return bricks[ib].cmatlist[iterm];
2119 void model::update_brick(
size_type ib, build_version version)
const {
2120 const brick_description &brick = bricks[ib];
2121 bool cplx =
is_complex() && brick.pbr->is_complex();
2122 bool tobecomputed = brick.terms_to_be_computed
2123 || brick.pbr->is_to_be_computed_each_time()
2124 || !(brick.pbr->is_linear());
2127 if (!tobecomputed ) {
2128 for (
size_type i = 0; i < brick.vlist.size() && !tobecomputed; ++i) {
2129 var_description &vd = variables[brick.vlist[i]];
2130 if (vd.v_num > brick.v_num) tobecomputed =
true;
2135 for (
size_type i = 0; i < brick.dlist.size() && !tobecomputed; ++i) {
2136 var_description &vd = variables[brick.dlist[i]];
2137 if (vd.v_num > brick.v_num || vd.v_num_data[vd.default_iter] > brick.v_num) {
2138 tobecomputed =
true;
2139 version = build_version(version | BUILD_ON_DATA_CHANGE);
2144 if (!tobecomputed ) {
2145 for (
size_type i = 0; i < brick.mims.size() && !tobecomputed; ++i) {
2146 if (brick.mims[i]->version_number() > brick.v_num) tobecomputed =
true;
2151 brick.external_load = scalar_type(0);
2153 if (!(brick.pdispatch))
2154 { brick_call(ib, version, 0); }
2157 brick.pdispatch->asm_complex_tangent_terms
2158 (*
this, ib, brick.cmatlist, brick.cveclist, brick.cveclist_sym,
2161 brick.pdispatch->asm_real_tangent_terms
2162 (*
this, ib, brick.rmatlist, brick.rveclist, brick.rveclist_sym,
2165 brick.v_num = act_counter();
2168 if (brick.pbr->is_linear()) brick.terms_to_be_computed =
false;
2175 const brick_description &brick = bricks[ib];
2176 if (brick.pbr->is_linear()) {
2178 bool cplx =
is_complex() && brick.pbr->is_complex();
2180 for (
size_type j = 0; j < brick.tlist.size(); ++j) {
2181 const term_description &term = brick.tlist[j];
2182 bool isg = term.is_global;
2185 size_type n_iter_1 = n_iter, n_iter_2 = n_iter;
2187 n_iter_1 = variables[term.var1].default_iter;
2188 if (term.is_matrix_term)
2189 n_iter_2 = variables[term.var2].default_iter;
2194 if (term.is_matrix_term) {
2197 model_complex_plain_vector V(nbgdof);
2198 for (VAR_SET::iterator it = variables.begin();
2199 it != variables.end(); ++it)
2200 if (it->second.is_variable) {
2202 ? it->second.default_iter : n_iter;
2203 gmm::copy(it->second.complex_value[n_iter_i],
2204 gmm::sub_vector(V, it->second.I));
2208 gmm::scaled(V, complex_type(-1)),
2209 brick.cveclist[ind_data][j]);
2213 gmm::scaled(variables[term.var2].complex_value[n_iter_2],
2215 brick.cveclist[ind_data][j]);
2219 model_real_plain_vector V(nbgdof);
2220 for (VAR_SET::iterator it = variables.begin();
2221 it != variables.end(); ++it)
2222 if (it->second.is_variable) {
2224 ? it->second.default_iter : n_iter;
2225 gmm::copy(it->second.real_value[n_iter_i],
2226 gmm::sub_vector(V, it->second.I));
2229 (brick.rmatlist[j], gmm::scaled(V, scalar_type(-1)),
2230 brick.rveclist[ind_data][j]);
2234 gmm::scaled(variables[term.var2].real_value[n_iter_2],
2235 scalar_type(-1)), brick.rveclist[ind_data][j]);
2238 if (term.is_symmetric && term.var1.compare(term.var2)) {
2241 (gmm::conjugated(brick.cmatlist[j]),
2242 gmm::scaled(variables[term.var1].complex_value[n_iter_1],
2244 brick.cveclist_sym[ind_data][j]);
2247 (gmm::transposed(brick.rmatlist[j]),
2248 gmm::scaled(variables[term.var1].real_value[n_iter_1],
2250 brick.rveclist_sym[ind_data][j]);
2257 void model::update_affine_dependent_variables() {
2258 for (VAR_SET::iterator it = variables.begin(); it != variables.end(); ++it)
2259 if (it->second.is_affine_dependent) {
2260 VAR_SET::iterator it2 = variables.find(it->second.org_name);
2261 if (it->second.size() != it2->second.size())
2262 it->second.set_size();
2263 if (it->second.is_complex) {
2264 gmm::add(gmm::scaled(it2->second.complex_value[0],
2265 complex_type(it->second.alpha)),
2266 it->second.affine_complex_value,
2267 it->second.complex_value[0]);
2269 gmm::add(gmm::scaled(it2->second.real_value[0], it->second.alpha),
2270 it->second.affine_real_value, it->second.real_value[0]);
2272 it->second.v_num = std::max(it->second.v_num, it2->second.v_num);
2273 for (
size_type i = 0; i < it->second.n_iter; ++i)
2275 it->second.v_num_data[i] = std::max(it->second.v_num_data[i],
2276 it2->second.v_num_data[i]);
2287 GMM_ASSERT1(mf,
"Works only with fem variables.");
2291 for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) {
2292 brick_description &brick = bricks[ib];
2294 bool detected =
false;
2295 for (
size_type i = 0; i < brick.vlist.size(); ++i)
2296 if (brick.vlist[i].compare(varname) == 0)
2297 { detected =
true;
break; }
2299 if (detected && brick.mims.size()) {
2301 for (
auto &pmim : brick.mims)
2304 pmim->linked_mesh()));
2305 GMM_ASSERT1(ifo >= 0,
2306 "The given region is only partially covered by "
2307 "region of brick \"" << brick.pbr->brick_name()
2308 <<
"\". Please subdivise the region");
2310 std::string expr = brick.pbr->declare_volume_assembly_string
2311 (*
this, ib, brick.vlist, brick.dlist);
2313 ga_workspace workspace(*
this);
2314 size_type order = workspace.add_expression
2315 (expr, dummy_mim, region);
2316 GMM_ASSERT1(order <= 1,
"Wrong order for a Neumann term");
2317 expr = workspace.extract_Neumann_term(varname);
2320 result +=
" + " + workspace.extract_Neumann_term(varname);
2322 result = workspace.extract_Neumann_term(varname);
2334 GMM_ASSERT1(version != BUILD_ON_DATA_CHANGE,
2335 "Invalid assembly version BUILD_ON_DATA_CHANGE");
2336 GMM_ASSERT1(version != BUILD_WITH_LIN,
2337 "Invalid assembly version BUILD_WITH_LIN");
2338 GMM_ASSERT1(version != BUILD_WITH_INTERNAL,
2339 "Invalid assembly version BUILD_WITH_INTERNAL");
2341 #if GETFEM_PARA_LEVEL > 1
2342 double t_ref = MPI_Wtime();
2344 MPI_Comm_rank(MPI_COMM_WORLD, &rk);
2345 MPI_Comm_size(MPI_COMM_WORLD, &nbp);
2348 context_check();
if (act_size_to_be_done) actualize_sizes();
2350 if (version & BUILD_MATRIX) gmm::clear(cTM);
2351 if (version & BUILD_RHS) gmm::clear(crhs);
2354 if (version & BUILD_MATRIX) gmm::clear(rTM);
2355 if (version & BUILD_RHS) gmm::clear(rrhs);
2357 clear_dof_constraints();
2358 generic_expressions.clear();
2359 update_affine_dependent_variables();
2361 if (version & BUILD_RHS) approx_external_load_ = scalar_type(0);
2363 for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) {
2365 brick_description &brick = bricks[ib];
2368 bool auto_disabled_brick =
true;
2369 for (
size_type j = 0; j < brick.vlist.size(); ++j) {
2371 auto_disabled_brick =
false;
2373 if (auto_disabled_brick)
continue;
2375 update_brick(ib, version);
2377 bool cplx =
is_complex() && brick.pbr->is_complex();
2379 scalar_type coeff0 = scalar_type(1);
2380 if (brick.pdispatch) coeff0 = brick.matrix_coeff;
2384 for (
size_type j = 0; j < brick.tlist.size(); ++j) {
2385 term_description &term = brick.tlist[j];
2386 bool isg = term.is_global, isprevious =
false;
2388 scalar_type alpha = coeff0, alpha1 = coeff0, alpha2 = coeff0;
2389 gmm::sub_interval I1(0,nbgdof), I2(0,nbgdof);
2390 var_description *var1=
nullptr, *var2=
nullptr;
2392 VAR_SET::iterator it1 = variables.find(term.var1);
2393 GMM_ASSERT1(it1 != variables.end(),
"Internal error");
2394 var1 = &(it1->second);
2395 GMM_ASSERT1(var1->is_variable,
"Assembly of data not allowed");
2397 if (term.is_matrix_term) {
2398 VAR_SET::iterator it2 = variables.find(term.var2);
2399 GMM_ASSERT1(it2 != variables.end(),
"Internal error");
2400 var2 = &(it2->second);
2402 if (!(var2->is_variable)) {
2403 std::string vorgname = sup_previous_and_dot_to_varname(term.var2);
2404 VAR_SET::iterator it3 = variables.find(vorgname);
2405 GMM_ASSERT1(it3->second.is_variable,
2406 "Assembly of data not allowed");
2410 alpha *= var1->alpha * var2->alpha;
2411 alpha1 *= var1->alpha;
2412 alpha2 *= var2->alpha;
2417 if (term.is_matrix_term && (version & BUILD_MATRIX) && !isprevious
2418 && (isg || (var1->is_enabled() && var2->is_enabled()))) {
2419 gmm::add(gmm::scaled(brick.cmatlist[j], alpha),
2420 gmm::sub_matrix(cTM, I1, I2));
2421 if (term.is_symmetric && I1.first() != I2.first())
2422 gmm::add(gmm::scaled(gmm::transposed(brick.cmatlist[j]), alpha),
2423 gmm::sub_matrix(cTM, I2, I1));
2425 if (version & BUILD_RHS) {
2427 if (isg || var1->is_enabled()) {
2428 if (brick.pdispatch)
2429 for (
size_type k = 0; k < brick.nbrhs; ++k)
2430 gmm::add(gmm::scaled(brick.cveclist[k][j],
2432 gmm::sub_vector(crhs, I1));
2434 gmm::add(gmm::scaled(brick.cveclist[0][j],
2435 complex_type(alpha1)),
2436 gmm::sub_vector(crhs, I1));
2438 if (term.is_matrix_term && brick.pbr->is_linear() &&
is_linear()) {
2439 if (var2->is_affine_dependent && var1->is_enabled())
2440 gmm::mult_add(brick.cmatlist[j],
2441 gmm::scaled(var2->affine_complex_value,
2442 complex_type(-alpha1)),
2443 gmm::sub_vector(crhs, I1));
2444 if (term.is_symmetric && I1.first() != I2.first()
2445 && var1->is_affine_dependent && var2->is_enabled())
2446 gmm::mult_add(gmm::conjugated(brick.cmatlist[j]),
2447 gmm::scaled(var1->affine_complex_value,
2448 complex_type(-alpha2)),
2449 gmm::sub_vector(crhs, I2));
2451 if (term.is_matrix_term && brick.pbr->is_linear()
2452 && (!
is_linear() || (version & BUILD_WITH_LIN))
2453 && var1->is_enabled())
2454 gmm::mult_add(brick.cmatlist[j],
2455 gmm::scaled(var2->complex_value[0],
2456 complex_type(-alpha1)),
2457 gmm::sub_vector(crhs, I1));
2458 if (term.is_symmetric && I1.first() != I2.first()
2459 && var2->is_enabled()) {
2460 if (brick.pdispatch)
2461 for (
size_type k = 0; k < brick.nbrhs; ++k)
2462 gmm::add(gmm::scaled(brick.cveclist_sym[k][j],
2464 gmm::sub_vector(crhs, I2));
2466 gmm::add(gmm::scaled(brick.cveclist_sym[0][j],
2467 complex_type(alpha2)),
2468 gmm::sub_vector(crhs, I2));
2469 if (brick.pbr->is_linear()
2470 && (!
is_linear() || (version & BUILD_WITH_LIN)))
2471 gmm::mult_add(gmm::conjugated(brick.cmatlist[j]),
2472 gmm::scaled(var1->complex_value[0],
2473 complex_type(-alpha2)),
2474 gmm::sub_vector(crhs, I2));
2478 if (term.is_matrix_term && (version & BUILD_MATRIX) && !isprevious
2479 && (isg || (var1->is_enabled() && var2->is_enabled()))) {
2480 gmm::add(gmm::scaled(brick.rmatlist[j], alpha),
2481 gmm::sub_matrix(cTM, I1, I2));
2482 if (term.is_symmetric && I1.first() != I2.first())
2483 gmm::add(gmm::scaled(gmm::transposed(brick.rmatlist[j]), alpha),
2484 gmm::sub_matrix(cTM, I2, I1));
2486 if (version & BUILD_RHS) {
2488 if (isg || var1->is_enabled()) {
2489 if (brick.pdispatch)
2490 for (
size_type k = 0; k < brick.nbrhs; ++k)
2491 gmm::add(gmm::scaled(brick.rveclist[k][j],
2493 gmm::sub_vector(crhs, I1));
2495 gmm::add(gmm::scaled(brick.rveclist[0][j], alpha1),
2496 gmm::sub_vector(crhs, I1));
2498 if (term.is_matrix_term && brick.pbr->is_linear() &&
is_linear()) {
2499 if (var2->is_affine_dependent && var1->is_enabled())
2500 gmm::mult_add(brick.rmatlist[j],
2501 gmm::scaled(var2->affine_complex_value,
2502 complex_type(-alpha1)),
2503 gmm::sub_vector(crhs, I1));
2504 if (term.is_symmetric && I1.first() != I2.first()
2505 && var1->is_affine_dependent && var2->is_enabled())
2506 gmm::mult_add(gmm::transposed(brick.rmatlist[j]),
2507 gmm::scaled(var1->affine_complex_value,
2508 complex_type(-alpha2)),
2509 gmm::sub_vector(crhs, I2));
2511 if (term.is_matrix_term && brick.pbr->is_linear()
2512 && (!
is_linear() || (version & BUILD_WITH_LIN))
2513 && var1->is_enabled())
2514 gmm::mult_add(brick.rmatlist[j],
2515 gmm::scaled(var2->complex_value[0],
2516 complex_type(-alpha1)),
2517 gmm::sub_vector(crhs, I1));
2518 if (term.is_symmetric && I1.first() != I2.first()
2519 && var2->is_enabled()) {
2520 if (brick.pdispatch)
2521 for (
size_type k = 0; k < brick.nbrhs; ++k)
2522 gmm::add(gmm::scaled(brick.rveclist_sym[k][j],
2524 gmm::sub_vector(crhs, I2));
2526 gmm::add(gmm::scaled(brick.rveclist_sym[0][j], alpha2),
2527 gmm::sub_vector(crhs, I2));
2529 if (brick.pbr->is_linear()
2530 && (!
is_linear() || (version & BUILD_WITH_LIN)))
2531 gmm::mult_add(gmm::transposed(brick.rmatlist[j]),
2532 gmm::scaled(var1->complex_value[0],
2533 complex_type(-alpha2)),
2534 gmm::sub_vector(crhs, I2));
2538 if (term.is_matrix_term && (version & BUILD_MATRIX) && !isprevious
2539 && (isg || (var1->is_enabled() && var2->is_enabled()))) {
2540 gmm::add(gmm::scaled(brick.rmatlist[j], alpha),
2541 gmm::sub_matrix(rTM, I1, I2));
2542 if (term.is_symmetric && I1.first() != I2.first())
2543 gmm::add(gmm::scaled(gmm::transposed(brick.rmatlist[j]), alpha),
2544 gmm::sub_matrix(rTM, I2, I1));
2546 if (version & BUILD_RHS) {
2548 auto vec_out1 = gmm::sub_vector(rrhs, I1);
2549 if (isg || var1->is_enabled()) {
2550 if (brick.pdispatch)
2551 for (
size_type k = 0; k < brick.nbrhs; ++k)
2552 gmm::add(gmm::scaled(brick.rveclist[k][j],
2556 gmm::add(gmm::scaled(brick.rveclist[0][j], alpha1),
2559 if (var1->is_enabled()
2560 && term.is_matrix_term && brick.pbr->is_linear()) {
2561 bool affine_contrib(
is_linear() && var2->is_affine_dependent);
2562 bool linear_contrib(!
is_linear() || (version & BUILD_WITH_LIN));
2563 const auto &matj = brick.rmatlist[j];
2564 const auto vec_affine2 = gmm::scaled(var2->affine_real_value,
2566 const auto vec_linear2 = gmm::scaled(var2->real_value[0],
2569 model_real_plain_vector vec_tmp1(I1.size(), 0.);
2571 gmm::mult(matj, vec_affine2, vec_tmp1);
2573 gmm::mult_add(matj, vec_linear2, vec_tmp1);
2574 MPI_SUM_VECTOR(vec_tmp1);
2575 gmm::add(vec_tmp1, vec_out1);
2578 gmm::mult_add(matj, vec_affine2, vec_out1);
2580 gmm::mult_add(matj, vec_linear2, vec_out1);
2585 if (term.is_symmetric && I1.first() != I2.first() &&
2586 var2->is_enabled()) {
2587 auto vec_out2 = gmm::sub_vector(rrhs, I2);
2588 if (brick.pdispatch)
2589 for (
size_type k = 0; k < brick.nbrhs; ++k)
2590 gmm::add(gmm::scaled(brick.rveclist_sym[k][j],
2594 gmm::add(gmm::scaled(brick.rveclist_sym[0][j], alpha2),
2596 if (term.is_matrix_term && brick.pbr->is_linear()) {
2597 bool affine_contrib(
is_linear() && var1->is_affine_dependent);
2598 bool linear_contrib(!
is_linear() || (version & BUILD_WITH_LIN));
2599 const auto matj_trans = gmm::transposed(brick.rmatlist[j]);
2600 const auto vec_affine1 = gmm::scaled(var1->affine_real_value,
2602 const auto vec_linear1 = gmm::scaled(var1->real_value[0],
2605 model_real_plain_vector vec_tmp2(I2.size(),0.);
2607 gmm::mult(matj_trans, vec_affine1, vec_tmp2);
2609 gmm::mult_add(matj_trans, vec_linear1, vec_tmp2);
2610 MPI_SUM_VECTOR(vec_tmp2);
2611 gmm::add(vec_tmp2, vec_out2);
2614 gmm::mult_add(matj_trans, vec_affine1, vec_out2);
2616 gmm::mult_add(matj_trans, vec_linear1, vec_out2);
2624 if (brick.pbr->is_linear())
2625 brick.terms_to_be_computed =
false;
2637 if (version & BUILD_RHS) approx_external_load_ += brick.external_load;
2640 if (version & BUILD_RHS && version & BUILD_WITH_INTERNAL) {
2643 gmm::sub_interval IP(0,gmm::vect_size(rrhs));
2644 gmm::fill(full_rrhs, 0.);
2645 gmm::copy(rrhs, gmm::sub_vector(full_rrhs, IP));
2649 if (generic_expressions.size()) {
2652 if (version & BUILD_RHS)
2653 GMM_TRACE2(
"Global generic assembly RHS");
2654 if (version & BUILD_MATRIX)
2655 GMM_TRACE2(
"Global generic assembly tangent term");
2658 auto add_assignments_and_expressions_to_workspace =
2659 [&](ga_workspace &workspace)
2661 for (
const auto &ad : assignments)
2662 workspace.add_assignment_expression
2663 (ad.varname, ad.expr, ad.region, ad.order, ad.before);
2664 for (
const auto &ge : generic_expressions)
2665 workspace.add_expression(ge.expr, ge.mim, ge.region,
2666 2, ge.secondary_domain);
2669 const bool with_internal = version & BUILD_WITH_INTERNAL
2671 model_real_sparse_matrix intern_mat;
2672 model_real_plain_vector res0,
2675 size_type full_size = gmm::vect_size(full_rrhs),
2676 primary_size = gmm::vect_size(rrhs);
2678 if ((version & BUILD_RHS) || (version & BUILD_MATRIX && with_internal))
2679 gmm::resize(res0, with_internal ? full_size : primary_size);
2680 if (version & BUILD_MATRIX && with_internal)
2681 gmm::resize(res1, full_size);
2683 if (version & BUILD_MATRIX) {
2684 if (with_internal) {
2685 gmm::resize(intern_mat, full_size, primary_size);
2686 gmm::resize(res1, full_size);
2692 if (version & BUILD_RHS) {
2695 ga_workspace workspace(*
this);
2696 add_assignments_and_expressions_to_workspace(workspace);
2697 workspace.set_assembled_vector(res0_distro);
2698 workspace.assembly(1, with_internal);
2699 if (with_internal) {
2700 gmm::copy(res0_distro.get(), res1_distro.get());
2701 gmm::add(gmm::scaled(full_rrhs, scalar_type(-1)),
2703 workspace.set_assembled_vector(res1_distro);
2704 workspace.set_internal_coupling_matrix(intern_mat_distro);
2706 workspace.set_assembled_matrix(tangent_matrix_distro);
2707 workspace.assembly(2, with_internal);
2712 ga_workspace workspace(*
this);
2713 add_assignments_and_expressions_to_workspace(workspace);
2714 if (with_internal) {
2715 gmm::copy(gmm::scaled(full_rrhs, scalar_type(-1)),
2717 workspace.set_assembled_vector(res1_distro);
2718 workspace.set_internal_coupling_matrix(intern_mat_distro);
2720 workspace.set_assembled_matrix(tangent_matrix_distro);
2721 workspace.assembly(2, with_internal);
2725 else if (version & BUILD_RHS) {
2728 ga_workspace workspace(*
this);
2729 add_assignments_and_expressions_to_workspace(workspace);
2730 workspace.set_assembled_vector(res0_distro);
2731 workspace.assembly(1, with_internal);
2735 if (version & BUILD_RHS) {
2736 gmm::scale(res0, scalar_type(-1));
2737 if (with_internal) {
2738 gmm::sub_interval IP(0,gmm::vect_size(rrhs));
2739 gmm::add(gmm::sub_vector(res0, IP), rrhs);
2740 gmm::add(res0, full_rrhs);
2742 gmm::add(res0, rrhs);
2745 if (version & BUILD_MATRIX && with_internal) {
2746 gmm::scale(res1, scalar_type(-1));
2747 gmm::sub_interval IP(0, primary_size),
2748 II(primary_size, full_size-primary_size);
2749 gmm::copy(gmm::sub_matrix(intern_mat, II, IP), internal_rTM);
2750 gmm::add(gmm::sub_vector(res1, IP), rrhs);
2751 gmm::copy(gmm::sub_vector(res1, II), internal_sol);
2756 if ((version & BUILD_RHS) || (version & BUILD_MATRIX)) {
2758 std::vector<size_type> dof_indices;
2759 std::vector<complex_type> dof_pr_values;
2760 std::vector<complex_type> dof_go_values;
2761 for (
const auto &keyval : complex_dof_constraints) {
2762 const gmm::sub_interval &I = interval_of_variable(keyval.first);
2764 for (
const auto &val : keyval.second) {
2765 dof_indices.push_back(val.first + I.first());
2766 dof_go_values.push_back(val.second);
2767 dof_pr_values.push_back(V[val.first]);
2771 if (dof_indices.size()) {
2772 gmm::sub_index SI(dof_indices);
2773 gmm::sub_interval II(0,
nb_dof());
2775 if (version & BUILD_RHS) {
2776 if (MPI_IS_MASTER())
2777 approx_external_load_ += gmm::vect_norm1(dof_go_values);
2779 if (is_symmetric_) {
2780 scalar_type valnorm = gmm::vect_norm2(dof_go_values);
2781 if (valnorm > scalar_type(0)) {
2782 GMM_ASSERT1(version & BUILD_MATRIX,
"Rhs only for a "
2783 "symmetric linear problem with dof "
2784 "constraint not allowed");
2785 model_complex_plain_vector vv(gmm::vect_size(crhs));
2786 gmm::mult(gmm::sub_matrix(cTM, II, SI), dof_go_values, vv);
2788 gmm::add(gmm::scaled(vv, scalar_type(-1)), crhs);
2791 gmm::copy(dof_go_values, gmm::sub_vector(crhs, SI));
2793 gmm::add(dof_go_values,
2794 gmm::scaled(dof_pr_values, complex_type(-1)),
2795 gmm::sub_vector(crhs, SI));
2798 if (version & BUILD_MATRIX) {
2799 gmm::clear(gmm::sub_matrix(cTM, SI, II));
2800 if (is_symmetric_) gmm::clear(gmm::sub_matrix(cTM, II, SI));
2802 if (MPI_IS_MASTER()) {
2803 for (
size_type i = 0; i < dof_indices.size(); ++i)
2804 cTM(dof_indices[i], dof_indices[i]) = complex_type(1);
2809 std::vector<size_type> dof_indices;
2810 std::vector<scalar_type> dof_pr_values;
2811 std::vector<scalar_type> dof_go_values;
2812 for (
const auto &keyval : real_dof_constraints) {
2813 const gmm::sub_interval &I = interval_of_variable(keyval.first);
2814 const model_real_plain_vector &V =
real_variable(keyval.first);
2815 for (
const auto &val : keyval.second) {
2816 dof_indices.push_back(val.first + I.first());
2817 dof_go_values.push_back(val.second);
2818 dof_pr_values.push_back(V[val.first]);
2822 #if GETFEM_PARA_LEVEL > 1
2823 GMM_ASSERT1(MPI_IS_MASTER() || (dof_indices.size() == 0),
2824 "Sorry, for the moment, the dof constraints have to be "
2825 "added on the master process only");
2826 size_type dof_indices_size = dof_indices.size();
2827 MPI_BCAST0_SCALAR(dof_indices_size);
2828 dof_indices.resize(dof_indices_size);
2829 MPI_BCAST0_VECTOR(dof_indices);
2830 dof_pr_values.resize(dof_indices_size);
2831 MPI_BCAST0_VECTOR(dof_pr_values);
2832 dof_go_values.resize(dof_indices_size);
2833 MPI_BCAST0_VECTOR(dof_go_values);
2836 if (dof_indices.size()) {
2837 gmm::sub_index SI(dof_indices);
2838 gmm::sub_interval II(0,
nb_dof());
2840 if (version & BUILD_RHS) {
2841 if (MPI_IS_MASTER())
2842 approx_external_load_ += gmm::vect_norm1(dof_go_values);
2844 if (is_symmetric_) {
2845 scalar_type valnorm = gmm::vect_norm2(dof_go_values);
2846 if (valnorm > scalar_type(0)) {
2847 GMM_ASSERT1(version & BUILD_MATRIX,
"Rhs only for a "
2848 "symmetric linear problem with dof "
2849 "constraint not allowed");
2850 model_real_plain_vector vv(gmm::vect_size(rrhs));
2851 gmm::mult(gmm::sub_matrix(rTM, II, SI), dof_go_values, vv);
2853 gmm::add(gmm::scaled(vv, scalar_type(-1)), rrhs);
2856 gmm::copy(dof_go_values, gmm::sub_vector(rrhs, SI));
2858 gmm::add(dof_go_values,
2859 gmm::scaled(dof_pr_values, scalar_type(-1)),
2860 gmm::sub_vector(rrhs, SI));
2863 if (version & BUILD_MATRIX) {
2864 gmm::clear(gmm::sub_matrix(rTM, SI, II));
2865 if (is_symmetric_) gmm::clear(gmm::sub_matrix(rTM, II, SI));
2867 if (MPI_IS_MASTER()) {
2868 for (
size_type i = 0; i < dof_indices.size(); ++i)
2869 rTM(dof_indices[i], dof_indices[i]) = scalar_type(1);
2876 if (version & BUILD_RHS) {
2879 MPI_BCAST0_SCALAR(approx_external_load_);
2882 #if GETFEM_PARA_LEVEL > 1
2884 if (MPI_IS_MASTER()) cout <<
"Assembly time " << MPI_Wtime()-t_ref << endl;
2893 return it->second.associated_mf();
2899 return it->second.passociated_mf();
2903 model::qdims_of_variable(
const std::string &name)
const {
2905 const mesh_fem *mf = it->second.passociated_mf();
2906 const im_data *imd = it->second.imd;
2909 bgeot::multi_index mi = mf->get_qdims();
2910 if (n > 1 || it->second.qdims.size() > 1) {
2912 if (mi.back() == 1) { mi.back() *= it->second.qdims[0]; ++i; }
2913 for (; i < it->second.qdims.size(); ++i)
2914 mi.push_back(it->second.qdims[i]);
2918 bgeot::multi_index mi = imd->tensor_size();
2919 if (n > 1 || it->second.qdims.size() > 1) {
2921 if (mi.back() == 1) { mi.back() *= it->second.qdims[0]; ++i; }
2922 for (; i < it->second.qdims.size(); ++i)
2923 mi.push_back(it->second.qdims[i]);
2927 return it->second.qdims;
2930 size_type model::qdim_of_variable(
const std::string &name)
const {
2932 const mesh_fem *mf = it->second.passociated_mf();
2933 const im_data *imd = it->second.imd;
2936 return mf->get_qdim() * n;
2938 return imd->tensor_size().total_size() * n;
2944 const gmm::sub_interval &
2945 model::interval_of_variable(
const std::string &name)
const {
2947 if (act_size_to_be_done) actualize_sizes();
2948 VAR_SET::const_iterator it = find_variable(name);
2949 return it->second.I;
2952 const model_real_plain_vector &
2958 const model_real_plain_vector &
2960 GMM_ASSERT1(!complex_version,
"This model is a complex one");
2961 GMM_ASSERT1(!
is_old(name),
"Please don't use Old_ prefix in combination "
2962 "with variable version");
2964 auto it = variables.find(name);
2965 GMM_ASSERT1(it != variables.end(),
"Undefined variable " << name);
2966 if (act_size_to_be_done && it->second.is_fem_dofs) {
2967 if (it->second.filter != VDESCRFILTER_NO)
2970 it->second.set_size();
2972 if (niter ==
size_type(-1)) niter = it->second.default_iter;
2973 GMM_ASSERT1(it->second.n_iter + it->second.n_temp_iter > niter,
2974 "Invalid iteration number " << niter <<
" for " << name);
2975 return it->second.real_value[niter];
2978 const model_complex_plain_vector &
2984 const model_complex_plain_vector &
2986 GMM_ASSERT1(complex_version,
"This model is a real one");
2987 GMM_ASSERT1(!
is_old(name),
"Please don't use Old_ prefix in combination with"
2988 " variable version");
2990 auto it = variables.find(name);
2991 GMM_ASSERT1(it!=variables.end(),
"Undefined variable " << name);
2992 if (act_size_to_be_done && it->second.is_fem_dofs) {
2993 if (it->second.filter != VDESCRFILTER_NO)
2996 it->second.set_size();
2998 if (niter ==
size_type(-1)) niter = it->second.default_iter;
2999 GMM_ASSERT1(it->second.n_iter + it->second.n_temp_iter > niter,
3000 "Invalid iteration number "
3001 << niter <<
" for " << name);
3002 return it->second.complex_value[niter];
3005 model_real_plain_vector &
3012 model_real_plain_vector &
3014 GMM_ASSERT1(!complex_version,
"This model is a complex one");
3015 GMM_ASSERT1(!
is_old(name),
"Please don't use Old_ prefix in combination with"
3016 " variable version");
3018 auto it = variables.find(name);
3019 GMM_ASSERT1(it!=variables.end(),
"Undefined variable " << name);
3020 if (act_size_to_be_done && it->second.is_fem_dofs) {
3021 if (it->second.filter != VDESCRFILTER_NO)
3024 it->second.set_size();
3026 if (niter ==
size_type(-1)) niter = it->second.default_iter;
3027 it->second.v_num_data[niter] = act_counter();
3028 GMM_ASSERT1(it->second.n_iter + it->second.n_temp_iter > niter,
3029 "Invalid iteration number "
3030 << niter <<
" for " << name);
3031 return it->second.real_value[niter];
3034 model_complex_plain_vector &
3040 model_complex_plain_vector &
3042 GMM_ASSERT1(complex_version,
"This model is a real one");
3043 GMM_ASSERT1(!
is_old(name),
"Please don't use Old_ prefix in combination with"
3044 " variable version");
3046 auto it = variables.find(name);
3047 GMM_ASSERT1(it!=variables.end(),
"Undefined variable " << name);
3048 if (act_size_to_be_done && it->second.is_fem_dofs) {
3049 if (it->second.filter != VDESCRFILTER_NO)
3052 it->second.set_size();
3054 if (niter ==
size_type(-1)) niter = it->second.default_iter;
3055 it->second.v_num_data[niter] = act_counter();
3056 GMM_ASSERT1(it->second.n_iter + it->second.n_temp_iter > niter,
3057 "Invalid iteration number "
3058 << niter <<
" for " << name);
3059 return it->second.complex_value[niter];
3062 model_real_plain_vector &
3063 model::set_real_constant_part(
const std::string &name)
const {
3064 GMM_ASSERT1(!complex_version,
"This model is a complex one");
3066 VAR_SET::iterator it = variables.find(name);
3067 GMM_ASSERT1(it != variables.end(),
"Undefined variable " << name);
3068 GMM_ASSERT1(it->second.is_affine_dependent,
3069 "Only for affine dependent variables");
3070 if (act_size_to_be_done && it->second.is_fem_dofs) {
3071 if (it->second.filter != VDESCRFILTER_NO)
3074 it->second.set_size();
3076 for (
auto &v_num : it->second.v_num_data) v_num = act_counter();
3077 return it->second.affine_real_value;
3080 model_complex_plain_vector &
3081 model::set_complex_constant_part(
const std::string &name)
const {
3082 GMM_ASSERT1(complex_version,
"This model is a real one");
3084 VAR_SET::iterator it = variables.find(name);
3085 GMM_ASSERT1(it!=variables.end(),
"Undefined variable " << name);
3086 if (act_size_to_be_done && it->second.is_fem_dofs) {
3087 if (it->second.filter != VDESCRFILTER_NO)
3090 it->second.set_size();
3092 for (
auto &v_num : it->second.v_num_data) v_num = act_counter();
3093 return it->second.affine_complex_value;
3098 const brick_description &brick = bricks[ind_brick];
3099 update_brick(ind_brick, model::BUILD_ALL);
3101 brick.pbr->check_stiffness_matrix_and_rhs(*
this, ind_brick, brick.tlist,
3102 brick.vlist, brick.dlist, brick.mims, brick.rmatlist,
3103 brick.rveclist[0], brick.rveclist_sym[0], brick.region);
3106 void model::clear() {
3108 active_bricks.clear();
3109 valid_bricks.clear();
3110 real_dof_constraints.clear();
3111 complex_dof_constraints.clear();
3113 rTM = model_real_sparse_matrix();
3114 cTM = model_complex_sparse_matrix();
3115 rrhs = model_real_plain_vector();
3116 crhs = model_complex_plain_vector();
3121 void virtual_brick::full_asm_real_tangent_terms_(
const model &md,
size_type ind_brick,
3122 const model::varnamelist &term_list,
3123 const model::varnamelist &var_list,
3124 const model::mimlist &mim_list,
3125 model::real_matlist &rmatlist,
3126 model::real_veclist &rveclist,
3127 model::real_veclist &rveclist_sym,
3128 size_type region, build_version version)
const
3131 rveclist, rveclist_sym, region, version);
3133 rveclist, rveclist_sym, region, version);
3135 rveclist, rveclist_sym, region, version);
3140 const model::termlist& tlist,
3141 const model::varnamelist &vl,
3142 const model::varnamelist &dl,
3143 const model::mimlist &mims,
3144 model::real_matlist &matl,
3145 model::real_veclist &rvc1,
3146 model::real_veclist &rvc2,
3148 const scalar_type TINY)
const
3150 std::cout<<
"******Verifying stiffnesses of *******"<<std::endl;
3151 std::cout<<
"*** "<<brick_name()<<std::endl;
3154 std::map<std::string,size_type> rhs_index;
3155 for(
size_type iterm=0;iterm<matl.size();iterm++)
3156 if (tlist[iterm].var1==tlist[iterm].var2) rhs_index[tlist[iterm].var1]=iterm;
3158 if (rhs_index.size()==0){
3159 GMM_WARNING0(
"*** cannot verify stiffness for this brick***");
3162 full_asm_real_tangent_terms_(md, s, vl, dl, mims, matl, rvc1, rvc2,
3163 rg, model::BUILD_MATRIX);
3164 for(
size_type iterm=0;iterm<matl.size();iterm++)
3167 std::cout<<std::endl;
3168 std::cout<<
" Stiffness["<<tlist[iterm].var1
3169 <<
","<<tlist[iterm].var2<<
"]:"<<std::endl;
3172 std::cout<<
" "<<tlist[iterm].var1<<
" has zero size. Skipping this term"<<std::endl;
3177 std::cout<<
" "<<tlist[iterm].var2<<
" has zero size. Skipping this term"<<std::endl;
3181 model_real_sparse_matrix SM(matl[iterm]);
3182 gmm::fill(rvc1[rhs_index[tlist[iterm].var1]], 0.0);
3183 full_asm_real_tangent_terms_(md, s, vl, dl, mims, matl, rvc1, rvc2,
3184 rg, model::BUILD_RHS);
3185 if (gmm::mat_euclidean_norm(matl[iterm])<1e-12){
3186 std::cout<<
" The assembled matrix is nearly zero, skipping."<<std::endl;
3189 model_real_plain_vector RHS0(rvc1[rhs_index[tlist[iterm].var1]]);
3192 model_real_sparse_matrix fdSM(matl[iterm].nrows(), matl[iterm].ncols());
3194 model_real_plain_vector& RHS1 = rvc1[rhs_index[tlist[iterm].var1]];
3196 scalar_type relative_tiny = gmm::vect_norminf(RHS1)*TINY;
3197 if (relative_tiny < TINY) relative_tiny = TINY;
3199 for (
size_type j = 0; j < matl[iterm].ncols(); j++){
3200 U[j] += relative_tiny;
3201 gmm::fill(RHS1, 0.0);
3202 full_asm_real_tangent_terms_(md, s, vl, dl, mims, matl, rvc1, rvc2,
3203 rg, model::BUILD_RHS);
3204 for (
size_type i = 0; i<matl[iterm].nrows(); i++)
3205 fdSM(i, j) = (RHS0[i] - RHS1[i]) / relative_tiny;
3206 U[j] -= relative_tiny;
3209 model_real_sparse_matrix diffSM(matl[iterm].nrows(),matl[iterm].ncols());
3210 gmm::add(SM,gmm::scaled(fdSM,-1.0),diffSM);
3211 scalar_type norm_error_euc
3212 = gmm::mat_euclidean_norm(diffSM)/gmm::mat_euclidean_norm(SM)*100;
3213 scalar_type norm_error_1
3214 = gmm::mat_norm1(diffSM)/gmm::mat_norm1(SM)*100;
3215 scalar_type norm_error_max
3216 = gmm::mat_maxnorm(diffSM)/gmm::mat_maxnorm(SM)*100;
3219 scalar_type nsym_norm_error_euc=0.0;
3220 scalar_type nsym_norm_error_1=0.0;
3221 scalar_type nsym_norm_error_max=0.0;
3222 if (tlist[iterm].var1==tlist[iterm].var2){
3223 model_real_sparse_matrix diffSMtransposed(matl[iterm].nrows(),matl[iterm].ncols());
3224 gmm::add(gmm::transposed(fdSM),gmm::scaled(fdSM,-1.0),diffSMtransposed);
3226 = gmm::mat_euclidean_norm(diffSMtransposed)/gmm::mat_euclidean_norm(fdSM)*100;
3228 = gmm::mat_norm1(diffSMtransposed)/gmm::mat_norm1(fdSM)*100;
3230 = gmm::mat_maxnorm(diffSMtransposed)/gmm::mat_maxnorm(fdSM)*100;
3234 if(rvc1[0].size()<8){
3235 std::cout <<
"RHS Stiffness Matrix: \n";
3236 std::cout <<
"------------------------\n";
3237 for(
size_type i=0; i < rvc1[iterm].size(); ++i){
3239 for(
size_type j=0; j < rvc1[iterm].size(); ++j){
3240 std::cout << fdSM(i,j) <<
" ";
3244 std::cout <<
"Analytical Stiffness Matrix: \n";
3245 std::cout <<
"------------------------\n";
3246 for(
size_type i=0; i < rvc1[iterm].size(); ++i){
3248 for(
size_type j=0; j < rvc1[iterm].size(); ++j){
3249 std::cout << matl[iterm](i,j) <<
" ";
3253 std::cout <<
"Vector U: \n";
3254 std::cout <<
"------------------------\n";
3255 for(
size_type i=0; i < rvc1[iterm].size(); ++i){
3262 <<
"\n\nfinite diff test error_norm_eucl: " << norm_error_euc <<
"%\n"
3263 <<
"finite diff test error_norm1: " << norm_error_1 <<
"%\n"
3264 <<
"finite diff test error_max_norm: " << norm_error_max <<
"%\n\n\n";
3266 if (tlist[iterm].var1==tlist[iterm].var2){
3268 <<
"Nonsymmetrical test error_norm_eucl: "<< nsym_norm_error_euc<<
"%\n"
3269 <<
"Nonsymmetrical test error_norm1: " << nsym_norm_error_1 <<
"%\n"
3270 <<
"Nonsymmetrical test error_max_norm: " << nsym_norm_error_max <<
"%"
3290 struct gen_source_term_assembly_brick :
public virtual_brick {
3292 std::string expr, directvarname, directdataname;
3293 model::varnamelist vl_test1;
3294 std::string secondary_domain;
3297 const model::varnamelist &,
3298 const model::varnamelist &,
3299 const model::mimlist &mims,
3300 model::real_matlist &,
3301 model::real_veclist &vecl,
3302 model::real_veclist &,
3304 build_version)
const override {
3305 GMM_ASSERT1(vecl.size() == vl_test1.size()
3306 + ((directdataname.size() == 0) ? 0 : 1),
"Wrong number "
3307 "of terms for Generic source term assembly brick ");
3308 GMM_ASSERT1(mims.size() == 1,
"Generic source term assembly brick "
3309 "needs one and only one mesh_im");
3310 GMM_TRACE2(
"Generic source term assembly");
3312 gmm::clear(vecl[0]);
3315 const mesh_im &mim = *mims[0];
3320 ga_workspace workspace(md, ga_workspace::inherit::ALL);
3321 GMM_TRACE2(name <<
": generic source term assembly");
3322 workspace.add_expression(expr, mim, rg, 1, secondary_domain);
3323 workspace.assembly(1);
3324 const auto &V=workspace.assembled_vector();
3325 for (
size_type i = 0; i < vl_test1.size(); ++i) {
3326 const auto &I=workspace.interval_of_variable(vl_test1[i]);
3327 gmm::copy(gmm::sub_vector(V, I), vecl[i]);
3331 if (directvarname.size()) {
3336 void real_post_assembly_in_serial(
const model &md,
size_type ib,
3337 const model::varnamelist &,
3338 const model::varnamelist &,
3339 const model::mimlist &,
3340 model::real_matlist &,
3341 model::real_veclist &vecl,
3342 model::real_veclist &,
3344 build_version)
const override {
3345 scalar_type el = scalar_type(0);
3346 for (
size_type i=0; i < vecl.size(); ++i) el += gmm::vect_norm1(vecl[i]);
3347 md.add_external_load(ib, el);
3350 virtual std::string declare_volume_assembly_string
3351 (
const model &,
size_type,
const model::varnamelist &,
3352 const model::varnamelist &)
const {
3353 return std::string();
3356 gen_source_term_assembly_brick(
const std::string &expr_,
3357 std::string brickname,
3358 const model::varnamelist &vl_test1_,
3359 const std::string &directvarname_,
3360 const std::string &directdataname_,
3361 const std::string &secdom)
3362 : vl_test1(vl_test1_), secondary_domain(secdom) {
3363 if (brickname.size() == 0)
3364 brickname =
"Generic source term assembly brick";
3366 set_flags(brickname,
true ,
3370 directvarname = directvarname_; directdataname = directdataname_;
3375 static bool check_compatibility_vl_test(model &md,
3376 const model::varnamelist vl_test) {
3377 model::varnamelist org;
3378 for (
size_type i = 0; i < vl_test.size(); ++i) {
3379 if (md.is_affine_dependent_variable(vl_test[i]))
3380 org.push_back(md.org_variable(vl_test[i]));
3382 for (
size_type i = 0; i < vl_test.size(); ++i)
3383 for (
size_type j = 0; j < org.size(); ++j)
3384 if (vl_test[i].compare(org[j]) == 0)
return false;
3389 (model &md,
const mesh_im &mim,
const std::string &expr,
size_type region,
3390 const std::string &brickname, std::string directvarname,
3391 const std::string &directdataname,
bool return_if_nonlin,
3392 const std::string &secondary_domain) {
3394 ga_workspace workspace(md);
3395 size_type order = workspace.add_expression(expr, mim, region, 1,
3397 GMM_ASSERT1(order <= 1,
"Wrong order for a source term");
3398 model::varnamelist vl, vl_test1, vl_test2, dl;
3399 bool is_lin = workspace.used_variables(vl, vl_test1, vl_test2, dl, 1);
3400 if (!is_lin && return_if_nonlin)
return size_type(-1);
3401 GMM_ASSERT1(is_lin,
"Nonlinear term");
3402 GMM_ASSERT1(check_compatibility_vl_test(md, vl_test1),
3403 "This brick do not support the assembly on both an affine "
3404 "dependent variable and its original variable. "
3405 "Split the brick.");
3407 if (directdataname.size()) {
3408 vl.push_back(directvarname);
3409 dl.push_back(directdataname);
3410 }
else directvarname =
"";
3412 pbrick pbr = std::make_shared<gen_source_term_assembly_brick>
3413 (expr, brickname, vl_test1, directvarname, directdataname,
3417 for (
size_type i = 0; i < vl_test1.size(); ++i)
3418 tl.push_back(model::term_description(vl_test1[i]));
3419 if (directdataname.size())
3420 tl.push_back(model::term_description(directvarname));
3422 return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
3427 const std::string &brickname,
const std::string &directvarname,
3428 const std::string &directdataname,
bool return_if_nonlin) {
3429 return add_source_term_(md, mim, expr, region, brickname, directvarname,
3430 directdataname, return_if_nonlin,
"");
3435 const std::string &secondary_domain,
3436 const std::string &brickname,
const std::string &directvarname,
3437 const std::string &directdataname,
bool return_if_nonlin) {
3438 return add_source_term_(md, mim, expr, region, brickname, directvarname,
3439 directdataname, return_if_nonlin, secondary_domain);
3448 struct gen_linear_assembly_brick :
public virtual_brick {
3452 model::varnamelist vl_test1, vl_test2;
3453 std::string secondary_domain;
3455 virtual void asm_real_tangent_terms(
const model &md,
size_type ib,
3456 const model::varnamelist &,
3457 const model::varnamelist &dl,
3458 const model::mimlist &mims,
3459 model::real_matlist &matl,
3460 model::real_veclist &,
3461 model::real_veclist &,
3463 build_version version)
const {
3464 GMM_ASSERT1(matl.size() == vl_test1.size(),
3465 "Wrong number of terms for Generic linear assembly brick");
3466 GMM_ASSERT1(mims.size() == 1,
3467 "Generic linear assembly brick needs one and only one "
3469 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0);
3470 for (
size_type i = 0; i < dl.size(); ++i)
3471 recompute_matrix = recompute_matrix ||
3472 md.is_var_newer_than_brick(dl[i], ib);
3474 if (recompute_matrix) {
3476 ga_workspace workspace(md, ga_workspace::inherit::ALL);
3477 workspace.add_expression(expr, *(mims[0]), region, 2, secondary_domain);
3478 GMM_TRACE2(name <<
": generic matrix assembly");
3479 workspace.assembly(2);
3480 const auto &R=workspace.assembled_matrix();
3481 for (
size_type i = 0; i < vl_test1.size(); ++i) {
3482 scalar_type alpha = scalar_type(1)
3483 / ( workspace.factor_of_variable(vl_test1[i]) *
3484 workspace.factor_of_variable(vl_test2[i]));
3485 const auto &I1=workspace.interval_of_variable(vl_test1[i]),
3486 &I2=workspace.interval_of_variable(vl_test2[i]);
3487 gmm::copy(gmm::scaled(gmm::sub_matrix(R, I1, I2), alpha),
3493 virtual std::string declare_volume_assembly_string
3494 (
const model &,
size_type,
const model::varnamelist &,
3495 const model::varnamelist &)
const {
3496 return is_lower_dim ? std::string() : expr;
3499 gen_linear_assembly_brick(
const std::string &expr_,
const mesh_im &mim,
3501 bool is_coer, std::string brickname,
3502 const model::varnamelist &vl_test1_,
3503 const model::varnamelist &vl_test2_,
3504 const std::string &secdom)
3505 : vl_test1(vl_test1_), vl_test2(vl_test2_), secondary_domain(secdom) {
3506 if (brickname.size() == 0) brickname =
"Generic linear assembly brick";
3508 is_lower_dim = mim.is_lower_dimensional();
3509 set_flags(brickname,
true ,
3516 static bool check_compatibility_vl_test(model &md,
3517 const model::varnamelist vl_test1,
3518 const model::varnamelist vl_test2) {
3519 for (
size_type i = 0; i < vl_test1.size(); ++i)
3520 for (
size_type j = 0; j < vl_test1.size(); ++j) {
3521 bool is1 = md.is_affine_dependent_variable(vl_test1[i]);
3522 bool is2 = md.is_affine_dependent_variable(vl_test2[i]);
3524 const std::string &org1
3525 = is1 ? md.org_variable(vl_test1[i]) : vl_test1[i];
3526 const std::string &org2
3527 = is2 ? md.org_variable(vl_test2[i]) : vl_test2[i];
3528 bool is1_bis = md.is_affine_dependent_variable(vl_test1[j]);
3529 bool is2_bis = md.is_affine_dependent_variable(vl_test2[j]);
3530 const std::string &org1_bis = is1_bis ? md.org_variable(vl_test1[j])
3532 const std::string &org2_bis = is2_bis ? md.org_variable(vl_test2[j])
3534 if (org1.compare(org1_bis) == 0 && org2.compare(org2_bis))
3544 (model &md,
const mesh_im &mim,
const std::string &expr,
size_type region,
3545 bool is_sym,
bool is_coercive,
const std::string &brickname,
3546 bool return_if_nonlin,
const std::string &secondary_domain) {
3548 ga_workspace workspace(md, ga_workspace::inherit::ALL);
3549 size_type order = workspace.add_expression(expr, mim, region,
3550 2, secondary_domain);
3551 model::varnamelist vl, vl_test1, vl_test2, dl;
3552 bool is_lin = workspace.used_variables(vl, vl_test1, vl_test2, dl, 2);
3554 if (!is_lin && return_if_nonlin)
return size_type(-1);
3555 GMM_ASSERT1(is_lin,
"Nonlinear term");
3556 if (order == 0) { is_coercive = is_sym =
true; }
3558 std::string const_expr= workspace.extract_constant_term(mim.linked_mesh());
3559 if (const_expr.size()) {
3561 (md, mim, const_expr, region, brickname+
" (source term)",
3562 "",
"",
false, secondary_domain);
3567 GMM_ASSERT1(check_compatibility_vl_test(md, vl_test1, vl_test2),
3568 "This brick do not support the assembly on both an affine "
3569 "dependent variable and its original variable. "
3570 "Split the brick.");
3572 if (vl_test1.size()) {
3573 pbrick pbr = std::make_shared<gen_linear_assembly_brick>
3574 (expr, mim, is_sym, is_coercive, brickname, vl_test1, vl_test2,
3577 for (
size_type i = 0; i < vl_test1.size(); ++i)
3578 tl.push_back(model::term_description(vl_test1[i], vl_test2[i],
false));
3580 return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
3587 bool is_sym,
bool is_coercive,
const std::string &brickname,
3588 bool return_if_nonlin) {
3589 return add_linear_term_(md, mim, expr, region, is_sym, is_coercive,
3590 brickname, return_if_nonlin,
"");
3595 const std::string &secondary_domain,
bool is_sym,
bool is_coercive,
3596 const std::string &brickname,
bool return_if_nonlin) {
3597 return add_linear_term_(md, mim, expr, region, is_sym, is_coercive,
3598 brickname, return_if_nonlin, secondary_domain);
3608 struct gen_nonlinear_assembly_brick :
public virtual_brick {
3612 std::string secondary_domain;
3614 virtual void real_post_assembly_in_serial(
const model &md,
size_type ,
3615 const model::varnamelist &,
3616 const model::varnamelist &,
3617 const model::mimlist &mims,
3618 model::real_matlist &,
3619 model::real_veclist &,
3620 model::real_veclist &,
3622 build_version)
const {
3623 GMM_ASSERT1(mims.size() == 1,
3624 "Generic linear assembly brick needs one and only one "
3626 md.add_generic_expression(expr, *(mims[0]), region, secondary_domain);
3629 virtual std::string declare_volume_assembly_string
3630 (
const model &,
size_type,
const model::varnamelist &,
3631 const model::varnamelist &)
const {
3636 gen_nonlinear_assembly_brick(
const std::string &expr_,
const mesh_im &mim,
3639 std::string brickname,
3640 const std::string &secdom) {
3641 if (brickname.size() == 0) brickname =
"Generic linear assembly brick";
3643 secondary_domain = secdom;
3644 is_lower_dim = mim.is_lower_dimensional();
3645 set_flags(brickname,
false ,
3653 (model &md,
const mesh_im &mim,
const std::string &expr,
size_type region,
3654 bool is_sym,
bool is_coercive,
const std::string &brickname,
3655 const std::string &secondary_domain) {
3657 ga_workspace workspace(md);
3658 size_type order = workspace.add_expression(expr, mim, region, 2,
3660 GMM_ASSERT1(order < 2,
"Order two test functions (Test2) are not allowed"
3661 " in assembly string for nonlinear terms");
3662 model::varnamelist vl, vl_test1, vl_test2, ddl, dl;
3663 workspace.used_variables(vl, vl_test1, vl_test2, ddl, order);
3664 for (
size_type i = 0; i < ddl.size(); ++i)
3665 if (md.is_true_data(ddl[i])) dl.push_back(ddl[i]);
3666 else vl.push_back(ddl[i]);
3667 if (order == 0) { is_coercive = is_sym =
true; }
3668 pbrick pbr = std::make_shared<gen_nonlinear_assembly_brick>
3669 (expr, mim, is_sym, is_coercive, brickname, secondary_domain);
3673 return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
3679 bool is_sym,
bool is_coercive,
const std::string &brickname) {
3680 return add_nonlinear_term_(md, mim, expr, region, is_sym, is_coercive,
3686 const std::string &secondary_domain,
bool is_sym,
bool is_coercive,
3687 const std::string &brickname) {
3688 return add_nonlinear_term_(md, mim, expr, region, is_sym, is_coercive,
3689 brickname, secondary_domain);
3700 struct generic_elliptic_brick :
public virtual_brick {
3702 virtual void asm_real_tangent_terms(
const model &md,
size_type ,
3703 const model::varnamelist &vl,
3704 const model::varnamelist &dl,
3705 const model::mimlist &mims,
3706 model::real_matlist &matl,
3707 model::real_veclist &,
3708 model::real_veclist &,
3710 build_version)
const {
3711 GMM_ASSERT1(matl.size() == 1,
3712 "Generic elliptic brick has one and only one term");
3713 GMM_ASSERT1(mims.size() == 1,
3714 "Generic elliptic brick need one and only one mesh_im");
3715 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
3716 "Wrong number of variables for generic elliptic brick");
3718 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
3719 const mesh &m = mf_u.linked_mesh();
3720 size_type N = m.dim(), Q = mf_u.get_qdim(), s = 1;
3721 const mesh_im &mim = *mims[0];
3722 const model_real_plain_vector *A = 0;
3723 const mesh_fem *mf_a = 0;
3724 mesh_region rg(region);
3725 m.intersect_with_mpi_region(rg);
3726 if (dl.size() > 0) {
3727 A = &(md.real_variable(dl[0]));
3728 mf_a = md.pmesh_fem_of_variable(dl[0]);
3729 s = gmm::vect_size(*A);
3730 if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
3734 GMM_TRACE2(
"Generic elliptic term assembly");
3739 (matl[0], mim, mf_u, *mf_a, *A, rg);
3742 (matl[0], mim, mf_u, *mf_a, *A, rg);
3747 (matl[0], mim, mf_u, rg);
3750 (matl[0], mim, mf_u, rg);
3751 if (A) gmm::scale(matl[0], (*A)[0]);
3753 }
else if (s == N*N) {
3757 (matl[0], mim, mf_u, *mf_a, *A, rg);
3760 (matl[0], mim, mf_u, *mf_a, *A, rg);
3764 (matl[0], mim, mf_u, *A, rg);
3767 (matl[0], mim, mf_u, *A, rg);
3769 }
else if (s == N*N*Q*Q) {
3772 (matl[0], mim, mf_u, *mf_a, *A, rg);
3775 (matl[0], mim, mf_u, *A, rg);
3777 GMM_ASSERT1(
false,
"Bad format generic elliptic brick coefficient");
3781 virtual void real_post_assembly_in_serial(
const model &md,
size_type,
3782 const model::varnamelist &,
3783 const model::varnamelist &dl,
3784 const model::mimlist &,
3785 model::real_matlist &,
3786 model::real_veclist &,
3787 model::real_veclist &,
3789 build_version)
const
3791 const model_real_plain_vector *A = 0;
3792 const mesh_fem *mf_a = 0;
3795 A = &(md.real_variable(dl[0]));
3796 mf_a = md.pmesh_fem_of_variable(dl[0]);
3800 virtual void complex_post_assembly_in_serial(
const model &md,
size_type,
3801 const model::varnamelist &,
3802 const model::varnamelist &dl,
3803 const model::mimlist &,
3804 model::complex_matlist &,
3805 model::complex_veclist &,
3806 model::complex_veclist &,
3808 build_version)
const
3810 const model_real_plain_vector *A = 0;
3811 const mesh_fem *mf_a = 0;
3814 A = &(md.real_variable(dl[0]));
3815 mf_a = md.pmesh_fem_of_variable(dl[0]);
3819 virtual scalar_type asm_complex_tangent_terms(
const model &md,
size_type,
3820 const model::varnamelist &vl,
3821 const model::varnamelist &,
3822 const model::mimlist &,
3823 model::complex_matlist &matl,
3824 model::complex_veclist &,
3825 model::complex_veclist &,
3827 const model_complex_plain_vector &U = md.complex_variable(vl[0]);
3828 return gmm::abs(gmm::vect_hp(matl[0], U, U)) / scalar_type(2);
3832 virtual void asm_complex_tangent_terms(
const model &md,
size_type,
3833 const model::varnamelist &vl,
3834 const model::varnamelist &dl,
3835 const model::mimlist &mims,
3836 model::complex_matlist &matl,
3837 model::complex_veclist &,
3838 model::complex_veclist &,
3840 build_version)
const {
3841 GMM_ASSERT1(matl.size() == 1,
3842 "Generic elliptic brick has one and only one term");
3843 GMM_ASSERT1(mims.size() == 1,
3844 "Generic elliptic brick need one and only one mesh_im");
3845 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
3846 "Wrong number of variables for generic elliptic brick");
3848 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
3849 const mesh &m = mf_u.linked_mesh();
3850 size_type N = m.dim(), Q = mf_u.get_qdim(), s = 1;
3851 const mesh_im &mim = *mims[0];
3852 const model_real_plain_vector *A = 0;
3853 const mesh_fem *mf_a = 0;
3854 mesh_region rg(region);
3855 m.intersect_with_mpi_region(rg);
3858 if (dl.size() > 0) {
3859 A = &(md.real_variable(dl[0]));
3860 mf_a = md.pmesh_fem_of_variable(dl[0]);
3861 s = gmm::vect_size(*A);
3862 if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
3866 GMM_TRACE2(
"Generic elliptic term assembly");
3871 (matl[0], mim, mf_u, *mf_a, *A, rg);
3874 (matl[0], mim, mf_u, *mf_a, *A, rg);
3879 (gmm::real_part(matl[0]), mim, mf_u, rg);
3882 (gmm::real_part(matl[0]), mim, mf_u, rg);
3883 if (A) gmm::scale(matl[0], (*A)[0]);
3885 }
else if (s == N*N) {
3889 (matl[0], mim, mf_u, *mf_a, *A, rg);
3892 (matl[0], mim, mf_u, *mf_a, *A, rg);
3896 (matl[0], mim, mf_u, *A, rg);
3899 (matl[0], mim, mf_u, *A, rg);
3901 }
else if (s == N*N*Q*Q) {
3904 (matl[0], mim, mf_u, *mf_a, *A, rg);
3907 (matl[0], mim, mf_u, *A, rg);
3910 "Bad format generic elliptic brick coefficient");
3913 generic_elliptic_brick() {
3914 set_flags(
"Generic elliptic",
true ,
3922 const std::string &varname,
3925 pbrick pbr = std::make_shared<generic_elliptic_brick>();
3927 tl.push_back(model::term_description(varname, varname,
true));
3928 return md.
add_brick(pbr, model::varnamelist(1, varname),
3929 model::varnamelist(), tl, model::mimlist(1, &mim),
3932 std::string test_varname
3933 =
"Test_" + sup_previous_and_dot_to_varname(varname);
3938 expr =
"Grad_"+varname+
".Grad_"+test_varname;
3940 expr =
"Grad_"+varname+
":Grad_"+test_varname;
3942 "Laplacian",
false);
3947 const std::string &varname,
3948 const std::string &dataname,
3951 pbrick pbr = std::make_shared<generic_elliptic_brick>();
3953 tl.push_back(model::term_description(varname, varname,
true));
3954 return md.
add_brick(pbr, model::varnamelist(1, varname),
3955 model::varnamelist(1, dataname), tl,
3956 model::mimlist(1, &mim), region);
3958 std::string test_varname
3959 =
"Test_" + sup_previous_and_dot_to_varname(varname);
3972 if (qdim_data != 1) {
3973 GMM_ASSERT1(qdim_data == gmm::sqr(dim),
3974 "Wrong data size for generic elliptic brick");
3975 expr =
"((Reshape("+dataname+
",meshdim,meshdim))*Grad_"+varname+
").Grad_"
3978 expr =
"(("+dataname+
")*Grad_"+varname+
").Grad_"+test_varname;
3981 if (qdim_data != 1) {
3982 if (qdim_data == gmm::sqr(dim))
3983 expr =
"((Reshape("+dataname+
",meshdim,meshdim))*Grad_"+varname+
"):Grad_"
3985 else if (qdim_data == gmm::sqr(gmm::sqr(dim)))
3986 expr =
"((Reshape("+dataname+
",meshdim,meshdim,meshdim,meshdim))*Grad_"
3987 +varname+
"):Grad_"+test_varname;
3988 else GMM_ASSERT1(
false,
"Wrong data size for generic elliptic brick");
3990 expr =
"(("+dataname+
")*Grad_"+varname+
"):Grad_"+test_varname;
3994 (md, mim, expr, region,
true,
true,
"Generic elliptic",
true);
3997 "Generic elliptic (nonlinear)");
4009 struct source_term_brick :
public virtual_brick {
4011 void asm_real_tangent_terms(
const model &md,
size_type ,
4012 const model::varnamelist &vl,
4013 const model::varnamelist &dl,
4014 const model::mimlist &mims,
4015 model::real_matlist &,
4016 model::real_veclist &vecl,
4017 model::real_veclist &,
4019 build_version)
const override {
4020 GMM_ASSERT1(vecl.size() == 1,
4021 "Source term brick has one and only one term");
4022 GMM_ASSERT1(mims.size() == 1,
4023 "Source term brick need one and only one mesh_im");
4024 GMM_ASSERT1(vl.size() == 1 && dl.size() > 0 && dl.size() <= 2,
4025 "Wrong number of variables for source term brick");
4027 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4028 const mesh_im &mim = *mims[0];
4029 const model_real_plain_vector &A = md.real_variable(dl[0]);
4030 const mesh_fem *mf_data = md.pmesh_fem_of_variable(dl[0]);
4031 mesh_region rg(region);
4032 mim.linked_mesh().intersect_with_mpi_region(rg);
4035 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4037 GMM_ASSERT1(mf_u.get_qdim() == s,
4038 dl[0] <<
": bad format of source term data. "
4039 "Detected dimension is " << s <<
" should be "
4042 GMM_TRACE2(
"Source term assembly");
4048 if (dl.size() > 1) gmm::add(md.real_variable(dl[1]), vecl[0]);
4051 void real_post_assembly_in_serial(
const model &md,
size_type ib,
4052 const model::varnamelist &,
4053 const model::varnamelist &,
4054 const model::mimlist &,
4055 model::real_matlist &,
4056 model::real_veclist &vecl,
4057 model::real_veclist &,
4059 build_version)
const override
4061 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4065 void asm_complex_tangent_terms(
const model &md,
size_type ,
4066 const model::varnamelist &vl,
4067 const model::varnamelist &dl,
4068 const model::mimlist &mims,
4069 model::complex_matlist &,
4070 model::complex_veclist &vecl,
4071 model::complex_veclist &,
4073 build_version)
const override {
4074 GMM_ASSERT1(vecl.size() == 1,
4075 "Source term brick has one and only one term");
4076 GMM_ASSERT1(mims.size() == 1,
4077 "Source term brick need one and only one mesh_im");
4078 GMM_ASSERT1(vl.size() == 1 && dl.size() > 0 && dl.size() <= 2,
4079 "Wrong number of variables for source term brick");
4081 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4082 const mesh_im &mim = *mims[0];
4083 const model_complex_plain_vector &A = md.complex_variable(dl[0]);
4084 const mesh_fem *mf_data = md.pmesh_fem_of_variable(dl[0]);
4085 mesh_region rg(region);
4086 mim.linked_mesh().intersect_with_mpi_region(rg);
4089 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4091 GMM_ASSERT1(mf_u.get_qdim() == s,
"Bad format of source term data");
4093 GMM_TRACE2(
"Source term assembly");
4099 if (dl.size() > 1) gmm::add(md.complex_variable(dl[1]), vecl[0]);
4102 void complex_post_assembly_in_serial(
const model &md,
4104 const model::varnamelist &,
4105 const model::varnamelist &,
4106 const model::mimlist &,
4107 model::complex_matlist &,
4108 model::complex_veclist &vecl,
4109 model::complex_veclist &,
4110 size_type, build_version)
const override
4112 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4117 source_term_brick() {
4118 set_flags(
"Source term",
true ,
4128 const std::string &varname,
4129 const std::string &dataexpr,
4131 const std::string &directdataname) {
4133 pbrick pbr = std::make_shared<source_term_brick>();
4135 tl.push_back(model::term_description(varname));
4136 model::varnamelist vdata(1, dataexpr);
4137 if (directdataname.size()) vdata.push_back(directdataname);
4138 return md.
add_brick(pbr, model::varnamelist(1, varname),
4139 vdata, tl, model::mimlist(1, &mim), region);
4141 std::string test_varname
4142 =
"Test_" + sup_previous_and_dot_to_varname(varname);
4147 expr =
"("+dataexpr+
")*"+test_varname;
4149 expr =
"("+dataexpr+
")."+test_varname;
4150 size_type ib = add_source_term_generic_assembly_brick
4151 (md, mim, expr, region,
"Source term", varname, directdataname,
true);
4154 "Source term (nonlinear)");
4155 if (directdataname.size())
4156 add_source_term_generic_assembly_brick
4157 (md, mim,
"", region,
"Source term", varname, directdataname);
4169 struct normal_source_term_brick :
public virtual_brick {
4171 void asm_real_tangent_terms(
const model &md,
size_type ,
4172 const model::varnamelist &vl,
4173 const model::varnamelist &dl,
4174 const model::mimlist &mims,
4175 model::real_matlist &,
4176 model::real_veclist &vecl,
4177 model::real_veclist &,
4179 build_version)
const override {
4180 GMM_ASSERT1(vecl.size() == 1,
4181 "Source term brick has one and only one term");
4182 GMM_ASSERT1(mims.size() == 1,
4183 "Source term brick need one and only one mesh_im");
4184 GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
4185 "Wrong number of variables for source term brick");
4187 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4188 const mesh_im &mim = *mims[0];
4189 const model_real_plain_vector &A = md.real_variable(dl[0]);
4190 const mesh_fem *mf_data = md.pmesh_fem_of_variable(dl[0]);
4191 mesh_region rg(region);
4192 mim.linked_mesh().intersect_with_mpi_region(rg);
4194 size_type s = gmm::vect_size(A), N = mf_u.linked_mesh().dim();
4195 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4197 GMM_ASSERT1(mf_u.get_qdim()*N == s,
4198 dl[0] <<
": bad format of normal source term data. "
4199 "Detected dimension is " << s <<
" should be "
4202 GMM_TRACE2(
"source term assembly");
4209 void real_post_assembly_in_serial(
const model &md,
size_type ib,
4210 const model::varnamelist &,
4211 const model::varnamelist &,
4212 const model::mimlist &,
4213 model::real_matlist &,
4214 model::real_veclist &vecl,
4215 model::real_veclist &,
4217 build_version)
const override {
4218 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4222 virtual void asm_complex_tangent_terms(
const model &md,
size_type ,
4223 const model::varnamelist &vl,
4224 const model::varnamelist &dl,
4225 const model::mimlist &mims,
4226 model::complex_matlist &,
4227 model::complex_veclist &vecl,
4228 model::complex_veclist &,
4230 build_version)
const {
4231 GMM_ASSERT1(vecl.size() == 1,
4232 "Source term brick has one and only one term");
4233 GMM_ASSERT1(mims.size() == 1,
4234 "Source term brick need one and only one mesh_im");
4235 GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
4236 "Wrong number of variables for source term brick");
4238 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4239 const mesh_im &mim = *mims[0];
4240 const model_complex_plain_vector &A = md.complex_variable(dl[0]);
4241 const mesh_fem *mf_data = md.pmesh_fem_of_variable(dl[0]);
4242 mesh_region rg(region);
4243 mim.linked_mesh().intersect_with_mpi_region(rg);
4245 size_type s = gmm::vect_size(A), N = mf_u.linked_mesh().dim();
4246 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4248 GMM_ASSERT1(s == mf_u.get_qdim()*N,
"Bad format of source term data");
4250 GMM_TRACE2(
"source term assembly");
4257 void complex_post_assembly_in_serial(
const model &md,
size_type ib,
4258 const model::varnamelist &,
4259 const model::varnamelist &,
4260 const model::mimlist &,
4261 model::complex_matlist &,
4262 model::complex_veclist &vecl,
4263 model::complex_veclist &,
4265 build_version)
const override {
4266 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4269 normal_source_term_brick() {
4270 set_flags(
"Normal source term",
true ,
4280 const std::string &varname,
4281 const std::string &dataexpr,
4284 pbrick pbr = std::make_shared<normal_source_term_brick>();
4286 tl.push_back(model::term_description(varname));
4287 model::varnamelist vdata(1, dataexpr);
4288 return md.
add_brick(pbr, model::varnamelist(1, varname),
4289 vdata, tl, model::mimlist(1, &mim), region);
4291 std::string test_varname
4292 =
"Test_" + sup_previous_and_dot_to_varname(varname);
4297 expr =
"(("+dataexpr+
").Normal)*"+test_varname;
4299 expr =
"(Reshape("+dataexpr+
",qdim("+varname
4300 +
"),meshdim)*Normal)."+test_varname;
4301 return add_source_term_generic_assembly_brick
4302 (md, mim, expr, region,
"Source term");
4315 struct Dirichlet_condition_brick :
public virtual_brick {
4318 bool normal_component;
4319 const mesh_fem *mf_mult_;
4325 virtual void asm_real_tangent_terms(
const model &md,
size_type ib,
4326 const model::varnamelist &vl,
4327 const model::varnamelist &dl,
4328 const model::mimlist &mims,
4329 model::real_matlist &matl,
4330 model::real_veclist &vecl,
4331 model::real_veclist &,
4333 build_version version)
const {
4334 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
4335 "Dirichlet condition brick has one and only one term");
4336 GMM_ASSERT1(mims.size() == 1,
4337 "Dirichlet condition brick need one and only one mesh_im");
4338 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() <= 3,
4339 "Wrong number of variables for Dirichlet condition brick");
4341 model_real_sparse_matrix& rB = rB_th;
4342 model_real_plain_vector& rV = rV_th;
4344 bool penalized = (vl.size() == 1);
4345 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4346 const mesh_fem &mf_mult = penalized ? (mf_mult_ ? *mf_mult_ : mf_u)
4347 : md.mesh_fem_of_variable(vl[1]);
4348 const mesh_im &mim = *mims[0];
4349 const model_real_plain_vector *A = 0, *COEFF = 0, *H = 0;
4350 const mesh_fem *mf_data = 0, *mf_H = 0;
4351 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
4352 || (penalized && md.is_var_newer_than_brick(dl[0], ib));
4355 COEFF = &(md.real_variable(dl[0]));
4356 GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
4357 "Data for coefficient should be a scalar");
4360 size_type s = 0, ind = (penalized ? 1 : 0);
4361 if (dl.size() > ind) {
4362 A = &(md.real_variable(dl[ind]));
4363 mf_data = md.pmesh_fem_of_variable(dl[ind]);
4364 s = gmm::vect_size(*A);
4365 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4366 size_type ss = ((normal_component) ? 1 : mf_u.get_qdim());
4367 GMM_ASSERT1(s == ss, dl[ind] <<
": bad format of "
4368 "Dirichlet data. Detected dimension is " << s
4369 <<
" should be " << ss);
4372 if (dl.size() > ind + 1) {
4373 GMM_ASSERT1(H_version,
4374 "Wrong number of data for Dirichlet condition brick");
4375 H = &(md.real_variable(dl[ind+1]));
4376 mf_H = md.pmesh_fem_of_variable(dl[ind+1]);
4377 s = gmm::vect_size(*A);
4379 s = s * mf_H->get_qdim() / mf_H->nb_dof();
4380 GMM_ASSERT1(mf_H->get_qdim() == 1,
"Implemented only for mf_H "
4381 "a scalar finite element method");
4383 GMM_ASSERT1(s = gmm::sqr(mf_u.get_qdim()),
4384 dl[ind+1] <<
": bad format of Dirichlet data. "
4385 "Detected dimension is " << s <<
" should be "
4386 <<
size_type(gmm::sqr(mf_u.get_qdim())));
4389 mesh_region rg(region);
4390 mim.linked_mesh().intersect_with_mpi_region(rg);
4392 if (recompute_matrix) {
4393 model_real_sparse_matrix *B = &(matl[0]);
4394 if (penalized && (&mf_mult != &mf_u)) {
4401 GMM_TRACE2(
"Mass term assembly for Dirichlet condition");
4402 if (H_version || normal_component) {
4403 ga_workspace workspace;
4404 gmm::sub_interval Imult(0, mf_mult.nb_dof()), Iu(0, mf_u.nb_dof());
4405 base_vector u(mf_u.nb_dof());
4406 base_vector mult(mf_mult.nb_dof());
4407 workspace.add_fem_variable(
"u", mf_u, Iu, u);
4408 workspace.add_fem_variable(
"mult", mf_mult, Imult, mult);
4409 auto expression = std::string{};
4411 if (mf_H) workspace.add_fem_constant(
"A", *mf_H, *H);
4412 else workspace.add_fixed_size_constant(
"A", *H);
4413 expression = (mf_u.get_qdim() == 1) ?
4414 "Test_mult . (A . Test2_u)"
4416 "Test_mult. (Reshape(A, qdim(u), qdim(u)) . Test2_u)";
4417 }
else if (normal_component) {
4418 expression =
"Test_mult . (Test2_u . Normal)";
4420 workspace.add_expression(expression, mim, rg);
4421 workspace.set_assembled_matrix(*B);
4422 workspace.assembly(2);
4427 if (penalized && (&mf_mult != &mf_u)) {
4428 GMM_ASSERT1(MPI_IS_MASTER(),
"Sorry, the penalized option "
4429 "filtered by a multiplier space is not parallelized");
4430 gmm::mult(gmm::transposed(rB), rB, matl[0]);
4431 gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
4432 }
else if (penalized) {
4433 gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
4437 if (dl.size() > ind) {
4438 GMM_TRACE2(
"Source term assembly for Dirichlet condition");
4440 if (penalized && (&mf_mult != &mf_u)) {
4454 if (penalized && (&mf_mult != &mf_u)) {
4455 gmm::mult(gmm::transposed(rB), rV, vecl[0]);
4456 gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
4457 rV = model_real_plain_vector();
4458 }
else if (penalized)
4459 gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
4464 void real_post_assembly_in_serial(
const model &md,
size_type ib,
4465 const model::varnamelist &,
4466 const model::varnamelist &,
4467 const model::mimlist &,
4468 model::real_matlist &,
4469 model::real_veclist &vecl,
4470 model::real_veclist &,
4472 build_version)
const override {
4473 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4476 virtual void asm_complex_tangent_terms(
const model &md,
size_type ib,
4477 const model::varnamelist &vl,
4478 const model::varnamelist &dl,
4479 const model::mimlist &mims,
4480 model::complex_matlist &matl,
4481 model::complex_veclist &vecl,
4482 model::complex_veclist &,
4484 build_version version)
const {
4485 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
4486 "Dirichlet condition brick has one and only one term");
4487 GMM_ASSERT1(mims.size() == 1,
4488 "Dirichlet condition brick need one and only one mesh_im");
4489 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() <= 3,
4490 "Wrong number of variables for Dirichlet condition brick");
4492 model_complex_sparse_matrix& cB = cB_th;
4493 model_complex_plain_vector& cV = cV_th;
4495 bool penalized = (vl.size() == 1);
4496 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4497 const mesh_fem &mf_mult = penalized ? (mf_mult_ ? *mf_mult_ : mf_u)
4498 : md.mesh_fem_of_variable(vl[1]);
4499 const mesh_im &mim = *mims[0];
4500 const model_complex_plain_vector *A = 0, *COEFF = 0, *H = 0;
4501 const mesh_fem *mf_data = 0, *mf_H = 0;
4502 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
4503 || (penalized && md.is_var_newer_than_brick(dl[0], ib));
4506 COEFF = &(md.complex_variable(dl[0]));
4507 GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
4508 "Data for coefficient should be a scalar");
4511 size_type s = 0, ind = (penalized ? 1 : 0);
4512 if (dl.size() > ind) {
4513 A = &(md.complex_variable(dl[ind]));
4514 mf_data = md.pmesh_fem_of_variable(dl[ind]);
4515 s = gmm::vect_size(*A);
4516 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4517 size_type ss = s * ((normal_component) ? mf_u.linked_mesh().dim() : 1);
4518 GMM_ASSERT1(mf_u.get_qdim() == ss,
4519 dl[ind] <<
": bad format of Dirichlet data. "
4520 "Detected dimension is " << ss <<
" should be "
4524 if (dl.size() > ind + 1) {
4525 GMM_ASSERT1(H_version,
4526 "Wrong number of data for Dirichlet condition brick");
4527 H = &(md.complex_variable(dl[ind+1]));
4528 mf_H = md.pmesh_fem_of_variable(dl[ind+1]);
4529 s = gmm::vect_size(*A);
4531 s = s * mf_H->get_qdim() / mf_H->nb_dof();
4532 GMM_ASSERT1(mf_H->get_qdim() == 1,
"Implemented only for mf_H "
4533 "a scalar finite element method");
4535 GMM_ASSERT1(s = gmm::sqr(mf_u.get_qdim()),
4536 dl[ind+1] <<
": bad format of Dirichlet data. "
4537 "Detected dimension is " << s <<
" should be "
4538 <<
size_type(gmm::sqr(mf_u.get_qdim())));
4541 mesh_region rg(region);
4542 mim.linked_mesh().intersect_with_mpi_region(rg);
4544 if (recompute_matrix) {
4545 model_complex_sparse_matrix *B = &(matl[0]);
4546 if (penalized && (&mf_mult != &mf_u)) {
4553 GMM_TRACE2(
"Mass term assembly for Dirichlet condition");
4555 if (mf_u.get_qdim() == 1)
4556 asm_real_or_complex_1_param_mat(*B, mim, mf_mult, mf_H, *H, rg,
4557 "(A*Test_u).Test2_u");
4559 asm_real_or_complex_1_param_mat(*B, mim, mf_mult, mf_H, *H, rg,
4560 "(Reshape(A,qdim(u),qdim(u))*Test2_u).Test_u");
4576 else if (normal_component) {
4577 ga_workspace workspace;
4578 gmm::sub_interval Imult(0, mf_mult.nb_dof()), Iu(0, mf_u.nb_dof());
4579 base_vector mult(mf_mult.nb_dof()), u(mf_u.nb_dof());
4580 workspace.add_fem_variable(
"mult", mf_mult, Imult, mult);
4581 workspace.add_fem_variable(
"u", mf_u, Iu, u);
4582 workspace.add_expression(
"Test_mult.(Test2_u.Normal)", mim, rg);
4583 model_real_sparse_matrix BB(mf_mult.nb_dof(), mf_u.nb_dof());
4584 workspace.set_assembled_matrix(BB);
4585 workspace.assembly(2);
4601 if (penalized && (&mf_mult != &mf_u)) {
4602 gmm::mult(gmm::transposed(cB), cB, matl[0]);
4603 gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
4604 }
else if (penalized) {
4605 gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
4609 if (dl.size() > ind) {
4610 GMM_TRACE2(
"Source term assembly for Dirichlet condition");
4612 if (penalized && (&mf_mult != &mf_u)) {
4626 if (penalized && (&mf_mult != &mf_u)) {
4627 gmm::mult(gmm::transposed(cB), cV, vecl[0]);
4628 gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
4629 cV = model_complex_plain_vector();
4630 }
else if (penalized)
4631 gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
4635 void complex_post_assembly_in_serial(
const model &md,
size_type ib,
4636 const model::varnamelist &,
4637 const model::varnamelist &,
4638 const model::mimlist &,
4639 model::complex_matlist &,
4640 model::complex_veclist &vecl,
4641 model::complex_veclist &,
4643 build_version)
const override {
4644 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4648 virtual std::string declare_volume_assembly_string
4649 (
const model &,
size_type,
const model::varnamelist &,
4650 const model::varnamelist &)
const {
4651 return std::string();
4654 Dirichlet_condition_brick(
bool penalized,
bool H_version_,
4655 bool normal_component_,
4656 const mesh_fem *mf_mult__ = 0) {
4657 mf_mult_ = mf_mult__;
4658 H_version = H_version_;
4659 normal_component = normal_component_;
4660 GMM_ASSERT1(!(H_version && normal_component),
"Bad Dirichlet version");
4661 set_flags(penalized ?
"Dirichlet with penalization brick"
4662 :
"Dirichlet with multipliers brick",
4672 const std::string &multname,
size_type region,
4673 const std::string &dataname) {
4674 pbrick pbr = std::make_shared<Dirichlet_condition_brick>(
false,
false,
false);
4676 tl.push_back(model::term_description(multname, varname,
true));
4677 model::varnamelist vl(1, varname);
4678 vl.push_back(multname);
4679 model::varnamelist dl;
4680 if (dataname.size()) dl.push_back(dataname);
4681 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4687 const std::string &dataname) {
4688 std::string multname = md.
new_name(
"mult_on_" + varname);
4691 (md, mim, varname, multname, region, dataname);
4697 const std::string &dataname) {
4702 (md, mim, varname, mf_mult, region, dataname);
4711 scalar_type penalisation_coeff,
size_type region,
4712 const std::string &dataname,
const mesh_fem *mf_mult) {
4713 std::string coeffname = md.
new_name(
"penalization_on_" + varname);
4719 pbrick pbr = std::make_shared<Dirichlet_condition_brick>
4720 (
true,
false,
false, mf_mult);
4722 tl.push_back(model::term_description(varname, varname,
true));
4723 model::varnamelist vl(1, varname);
4724 model::varnamelist dl(1, coeffname);
4725 if (dataname.size()) dl.push_back(dataname);
4726 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4731 const std::string &multname,
size_type region,
4732 const std::string &dataname) {
4733 pbrick pbr = std::make_shared<Dirichlet_condition_brick>(
false,
false,
true);
4735 tl.push_back(model::term_description(multname, varname,
true));
4736 model::varnamelist vl(1, varname);
4737 vl.push_back(multname);
4738 model::varnamelist dl;
4739 if (dataname.size()) dl.push_back(dataname);
4740 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4746 const std::string &dataname) {
4747 std::string multname = md.
new_name(
"mult_on_" + varname);
4750 (md, mim, varname, multname, region, dataname);
4756 const std::string &dataname) {
4760 (md, mim, varname, mf_mult, region, dataname);
4765 scalar_type penalisation_coeff,
size_type region,
4766 const std::string &dataname,
const mesh_fem *mf_mult) {
4767 std::string coeffname = md.
new_name(
"penalization_on_" + varname);
4773 pbrick pbr = std::make_shared<Dirichlet_condition_brick>
4774 (
true,
false,
true, mf_mult);
4776 tl.push_back(model::term_description(varname, varname,
true));
4777 model::varnamelist vl(1, varname);
4778 model::varnamelist dl(1, coeffname);
4779 if (dataname.size()) dl.push_back(dataname);
4780 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4786 const std::string &multname,
size_type region,
4787 const std::string &dataname,
const std::string &Hname) {
4788 pbrick pbr = std::make_shared<Dirichlet_condition_brick>(
false,
true,
false);
4790 tl.push_back(model::term_description(multname, varname,
true));
4791 model::varnamelist vl(1, varname);
4792 vl.push_back(multname);
4793 model::varnamelist dl;
4794 dl.push_back(dataname);
4795 dl.push_back(Hname);
4796 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4802 const std::string &dataname,
const std::string &Hname) {
4803 std::string multname = md.
new_name(
"mult_on_" + varname);
4806 (md, mim, varname, multname, region, dataname, Hname);
4812 const std::string &dataname,
const std::string &Hname) {
4817 (md, mim, varname, mf_mult, region, dataname, Hname);
4822 scalar_type penalisation_coeff,
size_type region,
4823 const std::string &dataname,
const std::string &Hname,
4825 std::string coeffname = md.
new_name(
"penalization_on_" + varname);
4831 pbrick pbr = std::make_shared<Dirichlet_condition_brick>
4832 (
true,
true,
false, mf_mult);
4834 tl.push_back(model::term_description(varname, varname,
true));
4835 model::varnamelist vl(1, varname);
4836 model::varnamelist dl(1, coeffname);
4837 dl.push_back(dataname);
4838 dl.push_back(Hname);
4839 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4843 scalar_type penalisation_coeff) {
4847 GMM_ASSERT1(gmm::vect_size(d)==1,
4848 "Wrong coefficient size, may be not a Dirichlet brick "
4849 "with penalization");
4850 d[0] = penalisation_coeff;
4854 GMM_ASSERT1(gmm::vect_size(d)==1,
4855 "Wrong coefficient size, may be not a Dirichlet brick "
4856 "with penalization");
4857 d[0] = penalisation_coeff;
4867 struct simplification_Dirichlet_condition_brick :
public virtual_brick {
4869 virtual void asm_real_tangent_terms(
const model &md,
size_type ,
4870 const model::varnamelist &vl,
4871 const model::varnamelist &dl,
4872 const model::mimlist &mims,
4873 model::real_matlist &matl,
4874 model::real_veclist &vecl,
4875 model::real_veclist &,
4877 build_version )
const {
4878 if (MPI_IS_MASTER()) {
4880 GMM_ASSERT1(vecl.size() == 0 && matl.size() == 0,
4881 "Dirichlet condition brick by simplification has no term");
4882 GMM_ASSERT1(mims.size() == 0,
4883 "Dirichlet condition brick by simplification need no "
4885 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
4886 "Wrong number of variables for Dirichlet condition brick "
4887 "by simplification");
4889 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4890 const model_real_plain_vector *A = 0;
4891 const mesh_fem *mf_data = 0;
4894 if (dl.size() == 1) {
4895 A = &(md.real_variable(dl[0]));
4896 mf_data = md.pmesh_fem_of_variable(dl[0]);
4899 GMM_ASSERT1(mf_data == &mf_u,
"Sorry, for this brick, the data has"
4900 " to be defined on the same f.e.m. as the unknown");
4902 s = gmm::vect_size(*A);
4903 GMM_ASSERT1(mf_u.get_qdim() == s,
": bad format of "
4904 "Dirichlet data. Detected dimension is " << s
4905 <<
" should be " <<
size_type(mf_u.get_qdim()));
4909 mesh_region rg(region);
4912 if (mf_u.get_qdim() > 1 || (!mf_data && A)) {
4913 for (mr_visitor i(rg, mf_u.linked_mesh()); !i.finished(); ++i) {
4914 pfem pf = mf_u.fem_of_element(i.cv());
4916 GMM_ASSERT1(pf->target_dim() == 1,
4917 "Intrinsically vectorial fems are not allowed");
4918 GMM_ASSERT1(mf_data || pf->is_lagrange(),
"Constant Dirichlet "
4919 "data allowed for lagrange fems only");
4924 dal::bit_vector dofs = mf_u.dof_on_region(rg);
4926 if (A && !mf_data) {
4927 GMM_ASSERT1(dofs.card() % s == 0,
"Problem with dof vectorization");
4930 for (dal::bv_visitor i(dofs); !i.finished(); ++i) {
4932 if (A) val = (mf_data ? (*A)[i] : (*A)[i%s]);
4933 md.add_real_dof_constraint(vl[0], i, val);
4938 virtual void asm_complex_tangent_terms(
const model &md,
size_type ,
4939 const model::varnamelist &vl,
4940 const model::varnamelist &dl,
4941 const model::mimlist &mims,
4942 model::complex_matlist &matl,
4943 model::complex_veclist &vecl,
4944 model::complex_veclist &,
4946 build_version )
const {
4947 if (MPI_IS_MASTER()) {
4948 GMM_ASSERT1(vecl.size() == 0 && matl.size() == 0,
4949 "Dirichlet condition brick by simplification has no term");
4950 GMM_ASSERT1(mims.size() == 0,
4951 "Dirichlet condition brick by simplification need no "
4953 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
4954 "Wrong number of variables for Dirichlet condition brick "
4955 "by simplification");
4957 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4958 const model_complex_plain_vector *A = 0;
4959 const mesh_fem *mf_data = 0;
4962 if (dl.size() == 1) {
4963 A = &(md.complex_variable(dl[0]));
4964 mf_data = md.pmesh_fem_of_variable(dl[0]);
4967 GMM_ASSERT1(mf_data == &mf_u,
"Sorry, for this brick, the data has"
4968 " to be define on the same f.e.m. than the unknown");
4970 s = gmm::vect_size(*A);
4971 GMM_ASSERT1(mf_u.get_qdim() == s,
": bad format of "
4972 "Dirichlet data. Detected dimension is " << s
4973 <<
" should be " <<
size_type(mf_u.get_qdim()));
4977 mesh_region rg(region);
4980 if (mf_u.get_qdim() > 1 || (!mf_data && A)) {
4981 for (mr_visitor i(rg, mf_u.linked_mesh()); !i.finished(); ++i) {
4982 pfem pf = mf_u.fem_of_element(i.cv());
4984 GMM_ASSERT1(pf->target_dim() == 1,
4985 "Intrinsically vectorial fems are not allowed");
4986 GMM_ASSERT1(mf_data || pf->is_lagrange(),
"Constant Dirichlet "
4987 "data allowed for lagrange fems only");
4992 dal::bit_vector dofs = mf_u.dof_on_region(rg);
4994 if (A && !mf_data) {
4995 GMM_ASSERT1(dofs.card() % s == 0,
"Problem with dof vectorization");
4998 for (dal::bv_visitor i(dofs); !i.finished(); ++i) {
4999 complex_type val(0);
5000 if (A) val = (mf_data ? (*A)[i] : (*A)[i%s]);
5001 md.add_complex_dof_constraint(vl[0], i, val);
5007 virtual std::string declare_volume_assembly_string
5008 (
const model &,
size_type,
const model::varnamelist &,
5009 const model::varnamelist &)
const {
5010 return std::string();
5013 simplification_Dirichlet_condition_brick() {
5014 set_flags(
"Dirichlet with simplification brick",
5024 size_type region,
const std::string &dataname) {
5025 pbrick pbr = std::make_shared<simplification_Dirichlet_condition_brick>();
5027 model::varnamelist vl(1, varname);
5028 model::varnamelist dl;
5029 if (dataname.size()) dl.push_back(dataname);
5030 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(), region);
5041 const std::string &Neumannterm,
5042 const std::string &datagamma0,
size_type region, scalar_type theta_,
5043 const std::string &datag) {
5044 std::string theta = std::to_string(theta_);
5046 ga_workspace workspace(md, ga_workspace::inherit::ALL);
5047 size_type order = workspace.add_expression(Neumannterm, mim, region, 1);
5048 GMM_ASSERT1(order == 0,
"Wrong expression of the Neumann term");
5049 bool is_lin = workspace.is_linear(1);
5051 std::string condition =
"("+varname + (datag.size() ?
"-("+datag+
"))":
")");
5052 std::string gamma =
"(("+datagamma0+
")*element_size)";
5053 std::string r =
"(1/"+gamma+
")";
5054 std::string expr =
"("+r+
"*"+condition+
"-("+Neumannterm+
")).Test_"+varname;
5055 if (theta_ != scalar_type(0)) {
5056 std::string derivative_Neumann = workspace.extract_order1_term(varname);
5057 if (derivative_Neumann.size())
5058 expr+=
"-"+theta+
"*"+condition+
".("+derivative_Neumann+
")";
5066 "Dirichlet condition with Nitsche's method");
5069 "Dirichlet condition with Nitsche's method");
5075 const std::string &Neumannterm,
5076 const std::string &datagamma0,
size_type region, scalar_type theta_,
5077 const std::string &datag) {
5078 std::string theta = std::to_string(theta_);
5080 ga_workspace workspace(md, ga_workspace::inherit::ALL);
5081 size_type order = workspace.add_expression(Neumannterm, mim, region, 1);
5082 GMM_ASSERT1(order == 0,
"Wrong expression of the Neumann term");
5083 bool is_lin = workspace.is_linear(1);
5085 std::string condition =
"("+varname+
".Normal"
5086 + (datag.size() ?
"-("+datag+
"))":
")");
5087 std::string gamma =
"(("+datagamma0+
")*element_size)";
5088 std::string r =
"(1/"+gamma+
")";
5089 std::string expr =
"("+r+
"*"+condition+
"-Normal.("+Neumannterm
5090 +
"))*(Normal.Test_"+varname+
")";
5091 if (theta_ != scalar_type(0)) {
5092 std::string derivative_Neumann = workspace.extract_order1_term(varname);
5093 if (derivative_Neumann.size())
5094 expr+=
"-"+theta+
"*"+condition+
"*Normal.("
5095 +derivative_Neumann+
")";
5099 "Dirichlet condition with Nitsche's method");
5102 "Dirichlet condition with Nitsche's method");
5108 const std::string &Neumannterm,
5109 const std::string &datagamma0,
size_type region, scalar_type theta_,
5110 const std::string &datag,
const std::string &dataH) {
5111 std::string theta = std::to_string(theta_);
5113 ga_workspace workspace(md, ga_workspace::inherit::ALL);
5114 size_type order = workspace.add_expression(Neumannterm, mim, region, 1);
5115 GMM_ASSERT1(order == 0,
"Wrong expression of the Neumann term");
5116 bool is_lin = workspace.is_linear(1);
5118 std::string condition =
"(("+dataH+
")*"+varname
5119 + (datag.size() ?
"-("+datag+
"))":
")");
5120 std::string gamma =
"(("+datagamma0+
")*element_size)";
5121 std::string r =
"(1/"+gamma+
")";
5122 std::string expr =
"("+r+
"*"+condition+
"-("+dataH+
")*("+Neumannterm
5123 +
"))*(("+dataH+
")*Test_"+varname+
")";
5124 if (theta_ != scalar_type(0)) {
5125 std::string derivative_Neumann = workspace.extract_order1_term(varname);
5126 if (derivative_Neumann.size())
5127 expr+=
"-"+theta+
"*"+condition+
"*(("+dataH+
")*("
5128 +derivative_Neumann+
"))";
5132 "Dirichlet condition with Nitsche's method");
5135 "Dirichlet condition with Nitsche's method");
5147 struct pointwise_constraints_brick :
public virtual_brick {
5149 mutable gmm::row_matrix<model_real_sparse_vector> rB;
5150 mutable gmm::row_matrix<model_complex_sparse_vector> cB;
5152 virtual void real_pre_assembly_in_serial(
const model &md,
size_type ib,
5153 const model::varnamelist &vl,
5154 const model::varnamelist &dl,
5155 const model::mimlist &mims,
5156 model::real_matlist &matl,
5157 model::real_veclist &vecl,
5158 model::real_veclist &,
5160 build_version version)
const {
5161 if (MPI_IS_MASTER()) {
5163 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5164 "Pointwize constraints brick has only one term");
5165 GMM_ASSERT1(mims.size() == 0,
5166 "Pointwize constraints brick does not need a mesh_im");
5167 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2,
5168 "Wrong number of variables for pointwize constraints brick");
5169 bool penalized = (vl.size() == 1);
5170 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5171 dim_type N = mf_u.linked_mesh().dim(), Q = mf_u.get_qdim(), ind_pt = 0;
5173 GMM_ASSERT1(dl.size() == dlsize || dl.size() == dlsize+1,
5174 "Wrong number of data for pointwize constraints brick");
5177 const model_real_plain_vector *COEFF = 0;
5179 COEFF = &(md.real_variable(dl[0]));
5181 GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
5182 "Data for coefficient should be a scalar");
5185 const model_real_plain_vector &PT = md.real_variable(dl[ind_pt]);
5186 size_type nb_co = gmm::vect_size(PT) / N;
5188 dim_type ind_unitv = dim_type((Q > 1) ? ind_pt+1 : 0);
5189 const model_real_plain_vector &unitv =md.real_variable(dl[ind_unitv]);
5190 GMM_ASSERT1((!ind_unitv || gmm::vect_size(unitv) == nb_co * Q),
5191 "Wrong size for vector of unit vectors");
5193 dim_type ind_rhs = dim_type((Q > 1) ? ind_pt+2 : ind_pt+1);
5194 if (dl.size() <
size_type(ind_rhs + 1)) ind_rhs = 0;
5195 const model_real_plain_vector &rhs = md.real_variable(dl[ind_rhs]);
5196 GMM_ASSERT1((!ind_rhs || gmm::vect_size(rhs) == nb_co),
5197 "Wrong size for vector of rhs");
5199 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
5200 || (penalized && (md.is_var_newer_than_brick(dl[ind_pt], ib)
5201 || md.is_var_newer_than_brick(dl[ind_unitv], ib)
5202 || md.is_var_newer_than_brick(dl[ind_rhs], ib)));
5204 if (recompute_matrix) {
5205 gmm::row_matrix<model_real_sparse_vector> BB(nb_co*Q, mf_u.nb_dof());
5208 dal::bit_vector dof_untouched;
5209 getfem::mesh_trans_inv mti(mf_u.linked_mesh());
5212 gmm::copy(gmm::sub_vector(PT, gmm::sub_interval(i*N, N)), pt);
5215 gmm::row_matrix<model_real_sparse_vector> &BBB = ((Q > 1) ? BB : rB);
5216 model_real_plain_vector vv;
5217 interpolation(mf_u, mti, vv, vv, BBB, 1, 1, &dof_untouched);
5218 GMM_ASSERT1(dof_untouched.card() == 0,
5219 "Pointwize constraints : some of the points are outside "
5220 "the mesh: " << dof_untouched);
5225 gmm::add(gmm::scaled(gmm::mat_row(BB, i*Q+q), unitv[i*Q+q]),
5226 gmm::mat_row(rB, i));
5229 gmm::mult(gmm::transposed(rB), rB, matl[0]);
5230 gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
5232 gmm::copy(rB, matl[0]);
5237 gmm::mult(gmm::transposed(rB), rhs, vecl[0]);
5238 gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
5240 else gmm::copy(rhs, vecl[0]);
5246 virtual void complex_pre_assembly_in_serial(
const model &md,
size_type ib,
5247 const model::varnamelist &vl,
5248 const model::varnamelist &dl,
5249 const model::mimlist &mims,
5250 model::complex_matlist &matl,
5251 model::complex_veclist &vecl,
5252 model::complex_veclist &,
5254 build_version version)
const {
5255 if (MPI_IS_MASTER()) {
5256 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5257 "Pointwize constraints brick only one term");
5258 GMM_ASSERT1(mims.size() == 0,
5259 "Pointwize constraints brick does not need a mesh_im");
5260 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2,
5261 "Wrong number of variables for pointwize constraints brick");
5262 bool penalized = (vl.size() == 1);
5263 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5264 dim_type N = mf_u.linked_mesh().dim(), Q = mf_u.get_qdim(), ind_pt = 0;
5266 GMM_ASSERT1(dl.size() == dlsize || dl.size() == dlsize+1,
5267 "Wrong number of data for pointwize constraints brick");
5270 const model_complex_plain_vector *COEFF = 0;
5272 COEFF = &(md.complex_variable(dl[0]));
5274 GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
5275 "Data for coefficient should be a scalar");
5278 const model_complex_plain_vector &PT = md.complex_variable(dl[ind_pt]);
5279 size_type nb_co = gmm::vect_size(PT) / N;
5281 dim_type ind_unitv = dim_type((Q > 1) ? ind_pt+1 : 0);
5282 const model_complex_plain_vector &unitv
5283 = md.complex_variable(dl[ind_unitv]);
5284 GMM_ASSERT1((!ind_unitv || gmm::vect_size(unitv) == nb_co * Q),
5285 "Wrong size for vector of unit vectors");
5287 dim_type ind_rhs = dim_type((Q > 1) ? ind_pt+2 : ind_pt+1);
5288 if (dl.size() <
size_type(ind_rhs + 1)) ind_rhs = 0;
5289 const model_complex_plain_vector &rhs
5290 = md.complex_variable(dl[ind_rhs]);
5291 GMM_ASSERT1((!ind_rhs || gmm::vect_size(rhs) == nb_co),
5292 "Wrong size for vector of rhs");
5294 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
5295 || (penalized && (md.is_var_newer_than_brick(dl[ind_pt], ib)
5296 || md.is_var_newer_than_brick(dl[ind_unitv], ib)
5297 || md.is_var_newer_than_brick(dl[ind_rhs], ib)));
5299 if (recompute_matrix) {
5300 gmm::row_matrix<model_complex_sparse_vector>
5301 BB(nb_co*Q,mf_u.nb_dof());
5303 dal::bit_vector dof_untouched;
5304 getfem::mesh_trans_inv mti(mf_u.linked_mesh());
5307 gmm::copy(gmm::real_part
5308 (gmm::sub_vector(PT, gmm::sub_interval(i*N, N))), pt);
5311 gmm::row_matrix<model_complex_sparse_vector> &BBB = ((Q > 1) ? BB :cB);
5312 model_complex_plain_vector vv;
5313 interpolation(mf_u, mti, vv, vv, BBB, 1, 1, &dof_untouched);
5314 GMM_ASSERT1(dof_untouched.card() == 0,
5315 "Pointwize constraints : some of the points are outside "
5316 "the mesh: " << dof_untouched);
5321 gmm::add(gmm::scaled(gmm::mat_row(BB, i*Q+q), unitv[i*Q+q]),
5322 gmm::mat_row(cB, i));
5326 gmm::mult(gmm::transposed(cB), cB, matl[0]);
5327 gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
5329 gmm::copy(cB, matl[0]);
5335 gmm::mult(gmm::transposed(cB), rhs, vecl[0]);
5336 gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
5338 else gmm::copy(rhs, vecl[0]);
5344 virtual std::string declare_volume_assembly_string
5345 (
const model &,
size_type,
const model::varnamelist &,
5346 const model::varnamelist &)
const {
5347 return std::string();
5350 pointwise_constraints_brick(
bool penalized) {
5351 set_flags(penalized ?
"Pointwise cosntraints with penalization brick"
5352 :
"Pointwise cosntraints with multipliers brick",
5363 scalar_type penalisation_coeff,
const std::string &dataname_pt,
5364 const std::string &dataname_unitv,
const std::string &dataname_val) {
5365 std::string coeffname = md.
new_name(
"penalization_on_" + varname);
5371 pbrick pbr = std::make_shared<pointwise_constraints_brick>(
true);
5373 tl.push_back(model::term_description(varname, varname,
true));
5374 model::varnamelist vl(1, varname);
5375 model::varnamelist dl(1, coeffname);
5376 dl.push_back(dataname_pt);
5378 if (mf_u.
get_qdim() > 1) dl.push_back(dataname_unitv);
5379 if (dataname_val.size() > 0) dl.push_back(dataname_val);
5385 const std::string &multname,
const std::string &dataname_pt,
5386 const std::string &dataname_unitv,
const std::string &dataname_val) {
5387 pbrick pbr = std::make_shared<pointwise_constraints_brick>(
false);
5389 tl.push_back(model::term_description(multname, varname,
true));
5390 model::varnamelist vl(1, varname);
5391 vl.push_back(multname);
5392 model::varnamelist dl(1, dataname_pt);
5394 if (mf_u.
get_qdim() > 1) dl.push_back(dataname_unitv);
5395 if (dataname_val.size() > 0) dl.push_back(dataname_val);
5400 (
model &md,
const std::string &varname,
const std::string &dataname_pt,
5401 const std::string &dataname_unitv,
const std::string &dataname_val) {
5402 std::string multname = md.
new_name(
"mult_on_" + varname);
5410 (md, varname, multname, dataname_pt, dataname_unitv, dataname_val);
5420 struct Helmholtz_brick :
public virtual_brick {
5422 virtual void asm_real_tangent_terms(
const model &md,
size_type,
5423 const model::varnamelist &vl,
5424 const model::varnamelist &dl,
5425 const model::mimlist &mims,
5426 model::real_matlist &matl,
5427 model::real_veclist &,
5428 model::real_veclist &,
5430 build_version)
const {
5431 GMM_ASSERT1(matl.size() == 1,
5432 "Helmholtz brick has one and only one term");
5433 GMM_ASSERT1(mims.size() == 1,
5434 "Helmholtz brick need one and only one mesh_im");
5435 GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
5436 "Wrong number of variables for Helmholtz brick");
5438 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5439 const mesh &m = mf_u.linked_mesh();
5441 GMM_ASSERT1(Q == 1,
"Helmholtz brick is only for scalar field, sorry.");
5442 const mesh_im &mim = *mims[0];
5443 const mesh_fem *mf_a = 0;
5444 mesh_region rg(region);
5445 m.intersect_with_mpi_region(rg);
5446 const model_real_plain_vector *A = &(md.real_variable(dl[0]));
5447 mf_a = md.pmesh_fem_of_variable(dl[0]);
5448 s = gmm::vect_size(*A);
5449 if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
5452 GMM_TRACE2(
"Stiffness matrix assembly for Helmholtz problem");
5453 gmm::clear(matl[0]);
5454 model_real_plain_vector A2(gmm::vect_size(*A));
5455 for (
size_type i=0; i < gmm::vect_size(*A); ++i)
5456 A2[i] = gmm::sqr((*A)[i]);
5462 GMM_ASSERT1(
false,
"Bad format Helmholtz brick coefficient");
5465 virtual void asm_complex_tangent_terms(
const model &md,
size_type,
5466 const model::varnamelist &vl,
5467 const model::varnamelist &dl,
5468 const model::mimlist &mims,
5469 model::complex_matlist &matl,
5470 model::complex_veclist &,
5471 model::complex_veclist &,
5473 build_version)
const {
5474 GMM_ASSERT1(matl.size() == 1,
5475 "Helmholtz brick has one and only one term");
5476 GMM_ASSERT1(mims.size() == 1,
5477 "Helmholtz brick need one and only one mesh_im");
5478 GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
5479 "Wrong number of variables for Helmholtz brick");
5481 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5482 const mesh &m = mf_u.linked_mesh();
5484 GMM_ASSERT1(Q == 1,
"Helmholtz brick is only for scalar field, sorry.");
5485 const mesh_im &mim = *mims[0];
5486 const mesh_fem *mf_a = 0;
5487 mesh_region rg(region);
5488 m.intersect_with_mpi_region(rg);
5489 const model_complex_plain_vector *A = &(md.complex_variable(dl[0]));
5490 mf_a = md.pmesh_fem_of_variable(dl[0]);
5491 s = gmm::vect_size(*A);
5492 if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
5495 GMM_TRACE2(
"Stiffness matrix assembly for Helmholtz problem");
5497 model_complex_plain_vector A2(gmm::vect_size(*A));
5498 for (
size_type i=0; i < gmm::vect_size(*A); ++i)
5499 A2[i] = gmm::sqr((*A)[i]);
5505 GMM_ASSERT1(
false,
"Bad format Helmholtz brick coefficient");
5509 set_flags(
"Helmholtz",
true ,
5517 const std::string &varname,
5518 const std::string &dataexpr,
5521 pbrick pbr = std::make_shared<Helmholtz_brick>();
5523 tl.push_back(model::term_description(varname, varname,
true));
5524 return md.
add_brick(pbr, model::varnamelist(1, varname),
5525 model::varnamelist(1, dataexpr), tl,
5526 model::mimlist(1, &mim), region);
5528 std::string test_varname
5529 =
"Test_" + sup_previous_and_dot_to_varname(varname);
5530 std::string expr =
"Grad_"+varname+
".Grad_"+test_varname
5531 +
" + sqr("+dataexpr+
")*"+varname+
"*"+test_varname;
5537 "Helmholtz (nonlinear)");
5550 struct Fourier_Robin_brick :
public virtual_brick {
5552 virtual void asm_real_tangent_terms(
const model &md,
size_type,
5553 const model::varnamelist &vl,
5554 const model::varnamelist &dl,
5555 const model::mimlist &mims,
5556 model::real_matlist &matl,
5557 model::real_veclist &,
5558 model::real_veclist &,
5560 build_version)
const {
5561 GMM_ASSERT1(matl.size() == 1,
5562 "Fourier-Robin brick has one and only one term");
5563 GMM_ASSERT1(mims.size() == 1,
5564 "Fourier-Robin brick need one and only one mesh_im");
5565 GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
5566 "Wrong number of variables for Fourier-Robin brick");
5568 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5569 const mesh &m = mf_u.linked_mesh();
5571 const mesh_im &mim = *mims[0];
5572 const mesh_fem *mf_a = 0;
5573 mesh_region rg(region);
5574 m.intersect_with_mpi_region(rg);
5575 const model_real_plain_vector *A = &(md.real_variable(dl[0]));
5576 mf_a = md.pmesh_fem_of_variable(dl[0]);
5577 s = gmm::vect_size(*A);
5578 if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
5579 GMM_ASSERT1(s == Q*Q,
5580 "Bad format Fourier-Robin brick coefficient");
5582 GMM_TRACE2(
"Fourier-Robin term assembly");
5583 gmm::clear(matl[0]);
5587 asm_homogeneous_qu_term(matl[0], mim, mf_u, *A, rg);
5590 virtual void asm_complex_tangent_terms(
const model &md,
size_type,
5591 const model::varnamelist &vl,
5592 const model::varnamelist &dl,
5593 const model::mimlist &mims,
5594 model::complex_matlist &matl,
5595 model::complex_veclist &,
5596 model::complex_veclist &,
5598 build_version)
const {
5599 GMM_ASSERT1(matl.size() == 1,
5600 "Fourier-Robin brick has one and only one term");
5601 GMM_ASSERT1(mims.size() == 1,
5602 "Fourier-Robin brick need one and only one mesh_im");
5603 GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
5604 "Wrong number of variables for Fourier-Robin brick");
5606 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5607 const mesh &m = mf_u.linked_mesh();
5609 const mesh_im &mim = *mims[0];
5610 const mesh_fem *mf_a = 0;
5611 mesh_region rg(region);
5612 m.intersect_with_mpi_region(rg);
5613 const model_complex_plain_vector *A = &(md.complex_variable(dl[0]));
5614 mf_a = md.pmesh_fem_of_variable(dl[0]);
5615 s = gmm::vect_size(*A);
5616 if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
5617 GMM_ASSERT1(s == Q*Q,
5618 "Bad format Fourier-Robin brick coefficient");
5620 GMM_TRACE2(
"Fourier-Robin term assembly");
5625 asm_homogeneous_qu_term(matl[0], mim, mf_u, *A, rg);
5628 Fourier_Robin_brick() {
5629 set_flags(
"Fourier Robin condition",
true ,
5638 const std::string &varname,
5639 const std::string &dataexpr,
5642 pbrick pbr = std::make_shared<Fourier_Robin_brick>();
5644 tl.push_back(model::term_description(varname, varname,
true));
5645 return md.
add_brick(pbr, model::varnamelist(1, varname),
5646 model::varnamelist(1, dataexpr), tl,
5647 model::mimlist(1, &mim), region);
5649 std::string test_varname
5650 =
"Test_" + sup_previous_and_dot_to_varname(varname);
5651 std::string expr =
"(("+dataexpr+
")*"+varname+
")."+test_varname;
5653 "Fourier-Robin",
true);
5656 "Fourier-Robin (nonlinear)");
5667 struct have_private_data_brick :
public virtual_brick {
5669 model_real_sparse_matrix rB;
5670 model_complex_sparse_matrix cB;
5671 model_real_plain_vector rL;
5672 model_complex_plain_vector cL;
5676 struct constraint_brick :
public have_private_data_brick {
5678 virtual void real_pre_assembly_in_serial(
const model &md,
size_type,
5679 const model::varnamelist &vl,
5680 const model::varnamelist &dl,
5681 const model::mimlist &mims,
5682 model::real_matlist &matl,
5683 model::real_veclist &vecl,
5684 model::real_veclist &,
5686 if (MPI_IS_MASTER()) {
5688 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5689 "Constraint brick has one and only one term");
5690 GMM_ASSERT1(mims.size() == 0,
5691 "Constraint brick need no mesh_im");
5692 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() <= 1,
5693 "Wrong number of variables for constraint brick");
5695 bool penalized = (vl.size() == 1);
5696 const model_real_plain_vector *COEFF = 0;
5698 bool has_data = (nameL.compare(
"") != 0);
5700 GMM_ASSERT1(nameL.compare(dl.back()) == 0 &&
5701 md.variable_exists(nameL) && md.is_data(nameL),
5703 const model_real_plain_vector &
5704 rrL = has_data ? md.real_variable(nameL) : rL;
5707 COEFF = &(md.real_variable(dl[0]));
5708 GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
5709 "Data for coefficient should be a scalar");
5711 gmm::mult(gmm::transposed(rB),
5712 gmm::scaled(rrL, gmm::abs((*COEFF)[0])), vecl[0]);
5713 gmm::mult(gmm::transposed(rB),
5714 gmm::scaled(rB, gmm::abs((*COEFF)[0])), matl[0]);
5716 gmm::copy(rrL, vecl[0]);
5717 gmm::copy(rB, matl[0]);
5722 virtual void complex_pre_assembly_in_serial(
const model &md,
size_type,
5723 const model::varnamelist &vl,
5724 const model::varnamelist &dl,
5725 const model::mimlist &mims,
5726 model::complex_matlist &matl,
5727 model::complex_veclist &vecl,
5728 model::complex_veclist &,
5730 build_version)
const {
5731 if (MPI_IS_MASTER()) {
5733 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5734 "Constraint brick has one and only one term");
5735 GMM_ASSERT1(mims.size() == 0,
5736 "Constraint brick need no mesh_im");
5737 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() <= 1,
5738 "Wrong number of variables for constraint brick");
5740 bool penalized = (vl.size() == 1);
5741 const model_complex_plain_vector *COEFF = 0;
5743 bool has_data = (nameL.compare(
"") != 0);
5745 GMM_ASSERT1(nameL.compare(dl.back()) == 0 &&
5746 md.variable_exists(nameL) && md.is_data(nameL),
5748 const model_complex_plain_vector &
5749 ccL = has_data ? md.complex_variable(nameL) : cL;
5752 COEFF = &(md.complex_variable(dl[0]));
5753 GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
5754 "Data for coefficient should be a scalar");
5756 gmm::mult(gmm::transposed(cB),
5757 gmm::scaled(ccL, gmm::abs((*COEFF)[0])), vecl[0]);
5758 gmm::mult(gmm::transposed(cB),
5759 gmm::scaled(cB, gmm::abs((*COEFF)[0])), matl[0]);
5761 gmm::copy(ccL, vecl[0]);
5762 gmm::copy(cB, matl[0]);
5767 virtual std::string declare_volume_assembly_string
5768 (
const model &,
size_type,
const model::varnamelist &,
5769 const model::varnamelist &)
const {
5770 return std::string();
5773 constraint_brick(
bool penalized) {
5774 set_flags(penalized ?
"Constraint with penalization brick"
5775 :
"Constraint with multipliers brick",
5784 model_real_sparse_matrix &set_private_data_brick_real_matrix
5786 pbrick pbr = md.brick_pointer(indbrick);
5787 md.touch_brick(indbrick);
5788 have_private_data_brick *p =
dynamic_cast<have_private_data_brick *
>
5789 (
const_cast<virtual_brick *
>(pbr.get()));
5790 GMM_ASSERT1(p,
"Wrong type of brick");
5794 model_real_plain_vector &set_private_data_brick_real_rhs
5796 pbrick pbr = md.brick_pointer(indbrick);
5797 md.touch_brick(indbrick);
5798 have_private_data_brick *p =
dynamic_cast<have_private_data_brick *
>
5799 (
const_cast<virtual_brick *
>(pbr.get()));
5800 GMM_ASSERT1(p,
"Wrong type of brick");
5801 if (p->nameL.compare(
"") != 0) GMM_WARNING1(
"Rhs already set by data name");
5805 model_complex_sparse_matrix &set_private_data_brick_complex_matrix
5807 pbrick pbr = md.brick_pointer(indbrick);
5808 md.touch_brick(indbrick);
5809 have_private_data_brick *p =
dynamic_cast<have_private_data_brick *
>
5810 (
const_cast<virtual_brick *
>(pbr.get()));
5811 GMM_ASSERT1(p,
"Wrong type of brick");
5815 model_complex_plain_vector &set_private_data_brick_complex_rhs
5817 pbrick pbr = md.brick_pointer(indbrick);
5818 md.touch_brick(indbrick);
5819 have_private_data_brick *p =
dynamic_cast<have_private_data_brick *
>
5820 (
const_cast<virtual_brick *
>(pbr.get()));
5821 GMM_ASSERT1(p,
"Wrong type of brick");
5822 if (p->nameL.compare(
"") != 0) GMM_WARNING1(
"Rhs already set by data name");
5826 void set_private_data_rhs
5827 (model &md,
size_type indbrick,
const std::string &varname) {
5828 pbrick pbr = md.brick_pointer(indbrick);
5829 md.touch_brick(indbrick);
5830 have_private_data_brick *p =
dynamic_cast<have_private_data_brick *
>
5831 (
const_cast<virtual_brick *
>(pbr.get()));
5832 GMM_ASSERT1(p,
"Wrong type of brick");
5833 if (p->nameL.compare(varname) != 0) {
5834 model::varnamelist dl = md.datanamelist_of_brick(indbrick);
5835 if (p->nameL.compare(
"") == 0) dl.push_back(varname);
5836 else dl.back() = varname;
5837 md.change_data_of_brick(indbrick, dl);
5842 size_type add_constraint_with_penalization
5843 (model &md,
const std::string &varname, scalar_type penalisation_coeff) {
5844 std::string coeffname = md.new_name(
"penalization_on_" + varname);
5845 md.add_fixed_size_data(coeffname, 1);
5846 if (md.is_complex())
5847 md.set_complex_variable(coeffname)[0] = penalisation_coeff;
5849 md.set_real_variable(coeffname)[0] = penalisation_coeff;
5850 pbrick pbr = std::make_shared<constraint_brick>(
true);
5852 tl.push_back(model::term_description(varname, varname,
true));
5853 model::varnamelist vl(1, varname);
5854 model::varnamelist dl(1, coeffname);
5855 return md.add_brick(pbr, vl, dl, tl, model::mimlist(),
size_type(-1));
5858 size_type add_constraint_with_multipliers
5859 (model &md,
const std::string &varname,
const std::string &multname) {
5860 pbrick pbr = std::make_shared<constraint_brick>(
false);
5862 tl.push_back(model::term_description(multname, varname,
true));
5863 model::varnamelist vl(1, varname);
5864 vl.push_back(multname);
5865 model::varnamelist dl;
5866 return md.add_brick(pbr, vl, dl, tl, model::mimlist(),
size_type(-1));
5876 struct explicit_matrix_brick :
public have_private_data_brick {
5878 virtual void real_pre_assembly_in_serial(
const model &,
size_type,
5879 const model::varnamelist &vl,
5880 const model::varnamelist &dl,
5881 const model::mimlist &mims,
5882 model::real_matlist &matl,
5883 model::real_veclist &vecl,
5884 model::real_veclist &,
5886 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5887 "Explicit matrix has one and only one term");
5888 GMM_ASSERT1(mims.size() == 0,
"Explicit matrix need no mesh_im");
5889 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() == 0,
5890 "Wrong number of variables for explicit matrix brick");
5891 GMM_ASSERT1(gmm::mat_ncols(rB) == gmm::mat_ncols(matl[0]) &&
5892 gmm::mat_nrows(rB) == gmm::mat_nrows(matl[0]),
5893 "Explicit matrix brick dimension mismatch ("<<
5894 gmm::mat_ncols(rB)<<
"x"<<gmm::mat_nrows(rB)<<
") != ("<<
5895 gmm::mat_ncols(matl[0])<<
"x"<<gmm::mat_nrows(matl[0])<<
")");
5896 gmm::copy(rB, matl[0]);
5899 virtual void complex_pre_assembly_in_serial(
const model &,
size_type,
5900 const model::varnamelist &vl,
5901 const model::varnamelist &dl,
5902 const model::mimlist &mims,
5903 model::complex_matlist &matl,
5904 model::complex_veclist &vecl,
5905 model::complex_veclist &,
5907 build_version)
const {
5908 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5909 "Explicit matrix has one and only one term");
5910 GMM_ASSERT1(mims.size() == 0,
"Explicit matrix need no mesh_im");
5911 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() == 0,
5912 "Wrong number of variables for explicit matrix brick");
5913 gmm::copy(cB, matl[0]);
5916 virtual std::string declare_volume_assembly_string
5917 (
const model &,
size_type,
const model::varnamelist &,
5918 const model::varnamelist &)
const {
5919 return std::string();
5922 explicit_matrix_brick(
bool symmetric_,
bool coercive_) {
5923 set_flags(
"Explicit matrix brick",
5925 symmetric_ , coercive_ ,
5932 (model &md,
const std::string &varname1,
const std::string &varname2,
5933 bool issymmetric,
bool iscoercive) {
5934 pbrick pbr = std::make_shared<explicit_matrix_brick>(issymmetric,
5937 tl.push_back(model::term_description(varname1, varname2, issymmetric));
5938 model::varnamelist vl(1, varname1);
5939 vl.push_back(varname2);
5940 model::varnamelist dl;
5941 return md.add_brick(pbr, vl, dl, tl, model::mimlist(),
size_type(-1));
5950 struct explicit_rhs_brick :
public have_private_data_brick {
5952 virtual void real_pre_assembly_in_serial(
const model &,
size_type,
5953 const model::varnamelist &vl,
5954 const model::varnamelist &dl,
5955 const model::mimlist &mims,
5956 model::real_matlist &matl,
5957 model::real_veclist &vecl,
5958 model::real_veclist &,
5960 if (MPI_IS_MASTER()) {
5961 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5962 "Explicit rhs has one and only one term");
5963 GMM_ASSERT1(mims.size() == 0,
"Explicit rhs need no mesh_im");
5964 GMM_ASSERT1(vl.size() == 1 && dl.size() == 0,
5965 "Wrong number of variables for explicit rhs brick");
5966 gmm::copy(rL, vecl[0]);
5970 virtual void complex_pre_assembly_in_serial(
const model &,
size_type,
5971 const model::varnamelist &vl,
5972 const model::varnamelist &dl,
5973 const model::mimlist &mims,
5974 model::complex_matlist &matl,
5975 model::complex_veclist &vecl,
5976 model::complex_veclist &,
5978 build_version)
const {
5979 if (MPI_IS_MASTER()) {
5980 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5981 "Explicit rhs has one and only one term");
5982 GMM_ASSERT1(mims.size() == 0,
"Explicit rhs need no mesh_im");
5983 GMM_ASSERT1(vl.size() == 1 && dl.size() == 0,
5984 "Wrong number of variables for explicit rhs brick");
5985 gmm::copy(cL, vecl[0]);
5990 virtual std::string declare_volume_assembly_string
5991 (
const model &,
size_type,
const model::varnamelist &,
5992 const model::varnamelist &)
const {
5993 return std::string();
5996 explicit_rhs_brick() {
5997 set_flags(
"Explicit rhs brick",
6007 (model &md,
const std::string &varname) {
6008 pbrick pbr = std::make_shared<explicit_rhs_brick>();
6010 tl.push_back(model::term_description(varname));
6011 model::varnamelist vl(1, varname);
6012 model::varnamelist dl;
6013 return md.add_brick(pbr, vl, dl, tl, model::mimlist(),
size_type(-1));
6023 struct iso_lin_elasticity_new_brick :
public virtual_brick {
6025 std::string expr, dataname3;
6027 void asm_real_tangent_terms(
const model &md,
size_type ib,
6028 const model::varnamelist &vl,
6029 const model::varnamelist &dl,
6030 const model::mimlist &mims,
6031 model::real_matlist &matl,
6032 model::real_veclist &vecl,
6033 model::real_veclist &,
6035 build_version version)
const override {
6036 GMM_ASSERT1(vl.size() == 1,
"Linearized isotropic elasticity brick "
6037 "has one and only one variable");
6038 GMM_ASSERT1(matl.size() == 1,
"Linearized isotropic elasticity brick "
6039 "has one and only one term");
6040 GMM_ASSERT1(mims.size() == 1,
"Linearized isotropic elasticity brick "
6041 "needs one and only one mesh_im");
6043 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0);
6044 for (
size_type i = 0; i < dl.size(); ++i) {
6045 recompute_matrix = recompute_matrix ||
6046 md.is_var_newer_than_brick(dl[i], ib);
6049 if (recompute_matrix) {
6051 ga_workspace workspace(md, ga_workspace::inherit::ALL);
6052 workspace.add_expression(expr, *(mims[0]), region);
6053 GMM_TRACE2(name <<
": generic matrix assembly");
6054 workspace.assembly(2);
6055 scalar_type
alpha = scalar_type(1)
6056 / (workspace.factor_of_variable(vl[0]));
6057 const auto &R=workspace.assembled_matrix();
6058 gmm::sub_interval I = workspace.interval_of_variable(vl[0]);
6059 gmm::copy(gmm::scaled(gmm::sub_matrix(R, I, I), alpha),
6063 if (dataname3.size()) {
6069 gmm::scaled(md.real_variable(dataname3), scalar_type(-1)),
6075 void real_post_assembly_in_serial(
const model &md,
size_type ib,
6076 const model::varnamelist &,
6077 const model::varnamelist &,
6078 const model::mimlist &,
6079 model::real_matlist &,
6080 model::real_veclist &vecl,
6081 model::real_veclist &,
6083 build_version)
const override {
6084 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
6088 virtual std::string declare_volume_assembly_string
6089 (
const model &,
size_type,
const model::varnamelist &,
6090 const model::varnamelist &)
const {
6094 iso_lin_elasticity_new_brick(
const std::string &expr_,
6095 const std::string &dataname3_) {
6096 expr = expr_; dataname3 = dataname3_;
6097 set_flags(
"Linearized isotropic elasticity",
true ,
6107 const std::string &dataexpr1,
const std::string &dataexpr2,
6108 size_type region,
const std::string &dataname3) {
6109 std::string test_varname
6110 =
"Test_" + sup_previous_and_dot_to_varname(varname);
6112 std::string expr1 =
"((("+dataexpr1+
")*(Div_"+varname+
"-Div_"+dataname3
6113 +
"))*Id(meshdim)+(2*("+dataexpr2+
"))*(Sym(Grad_"+varname
6114 +
")-Sym(Grad_"+dataname3+
"))):Grad_" +test_varname;
6115 std::string expr2 =
"(Div_"+varname+
"*(("+dataexpr1+
")*Id(meshdim))"
6116 +
"+(2*("+dataexpr2+
"))*Sym(Grad_"+varname+
")):Grad_"+test_varname;
6119 model::varnamelist vl, dl;
6121 ga_workspace workspace(md, ga_workspace::inherit::ALL);
6122 workspace.add_expression(expr2, mim, region);
6123 model::varnamelist vl_test1, vl_test2;
6124 is_lin = workspace.used_variables(vl, vl_test1, vl_test2, dl, 2);
6127 pbrick pbr = std::make_shared<iso_lin_elasticity_new_brick>
6130 tl.push_back(model::term_description(varname,
6131 sup_previous_and_dot_to_varname(varname),
true));
6132 if (dataname3.size()) dl.push_back(dataname3);
6133 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
6135 return add_nonlinear_generic_assembly_brick
6136 (md, mim, dataname3.size() ? expr1 : expr2, region,
false,
false,
6137 "Linearized isotropic elasticity (with nonlinear dependance)");
6143 const std::string &data_E,
const std::string &data_nu,
6145 std::string test_varname
6146 =
"Test_" + sup_previous_and_dot_to_varname(varname);
6148 std::string mu =
"(("+data_E+
")/(2*(1+("+data_nu+
"))))";
6149 std::string lambda =
"(("+data_E+
")*("+data_nu+
")/((1+("+data_nu+
"))*(1-2*("
6151 std::string expr = lambda+
"*Div_"+varname+
"*Div_"+test_varname
6152 +
"+"+mu+
"*(Grad_"+varname+
"+Grad_"+varname+
"'):Grad_"+test_varname;
6156 ga_workspace workspace(md, ga_workspace::inherit::ALL);
6157 workspace.add_expression(expr, mim, region);
6158 is_lin = workspace.is_linear(2);
6162 "Linearized isotropic elasticity");
6165 (md, mim, expr, region,
false,
false,
6166 "Linearized isotropic elasticity (with nonlinear dependance)");
6172 const std::string &data_E,
const std::string &data_nu,
6174 std::string test_varname
6175 =
"Test_" + sup_previous_and_dot_to_varname(varname);
6178 GMM_ASSERT1(mfu,
"The variable should be a fem variable");
6181 std::string mu =
"(("+data_E+
")/(2*(1+("+data_nu+
"))))";
6182 std::string lambda =
"(("+data_E+
")*("+data_nu+
")/((1+("+data_nu+
"))*(1-2*("
6185 lambda =
"(("+data_E+
")*("+data_nu+
")/((1-sqr("+data_nu+
"))))";
6186 std::string expr = lambda+
"*Div_"+varname+
"*Div_"+test_varname
6187 +
"+"+mu+
"*(Grad_"+varname+
"+Grad_"+varname+
"'):Grad_"+test_varname;
6191 ga_workspace workspace(md, ga_workspace::inherit::ALL);
6192 workspace.add_expression(expr, mim, region);
6193 is_lin = workspace.is_linear(2);
6197 "Linearized isotropic elasticity");
6200 (md, mim, expr, region,
false,
false,
6201 "Linearized isotropic elasticity (with nonlinear dependance)");
6206 void compute_isotropic_linearized_Von_Mises_or_Tresca
6207 (model &md,
const std::string &varname,
const std::string &data_lambda,
6208 const std::string &data_mu,
const mesh_fem &mf_vm,
6209 model_real_plain_vector &VM,
bool tresca) {
6212 const mesh_fem &mf_u = md.mesh_fem_of_variable(varname);
6213 const mesh_fem *mf_lambda = md.pmesh_fem_of_variable(data_lambda);
6214 const model_real_plain_vector *lambda=&(md.real_variable(data_lambda));
6215 const mesh_fem *mf_mu = md.pmesh_fem_of_variable(data_mu);
6216 const model_real_plain_vector *mu = &(md.real_variable(data_mu));
6219 if (mf_lambda) sl = sl * mf_lambda->get_qdim() / mf_lambda->nb_dof();
6221 if (mf_mu) sm = sm * mf_mu->get_qdim() / mf_mu->nb_dof();
6223 GMM_ASSERT1(sl == 1 && sm == 1,
"Bad format for Lame coefficients");
6224 GMM_ASSERT1(mf_lambda == mf_mu,
6225 "The two Lame coefficients should be described on the same "
6226 "finite element method.");
6230 md.real_variable(varname), VM,
6231 *mf_lambda, *lambda,
6236 model_real_plain_vector LAMBDA(mf_lambda->nb_dof(), (*lambda)[0]);
6237 model_real_plain_vector MU(mf_lambda->nb_dof(), (*mu)[0]);
6239 md.real_variable(varname), VM,
6248 std::string sigma_d=
"("+data_mu+
")*(Grad_"+varname+
"+Grad_"+varname+
"')";
6249 std::string expr =
"sqrt(3/2)*Norm(Deviator("+sigma_d+
"))";
6250 ga_interpolation_Lagrange_fem(md, expr, mf_vm, VM);
6256 (
model &md,
const std::string &varname,
const std::string &data_E,
6257 const std::string &data_nu,
const mesh_fem &mf_vm,
6258 model_real_plain_vector &VM) {
6262 std::string mu =
"(("+data_E+
")/(2*(1+("+data_nu+
"))))";
6265 std::string sigma_d = mu+
"*(Grad_"+varname+
"+Grad_"+varname+
"')";
6266 std::string expr =
"sqrt(3/2)*Norm(Deviator("+sigma_d+
"))";
6267 ga_interpolation_Lagrange_fem(md, expr, mf_vm, VM);
6271 (
model &md,
const std::string &varname,
const std::string &data_E,
6272 const std::string &data_nu,
const mesh_fem &mf_vm,
6273 model_real_plain_vector &VM) {
6282 std::string mu =
"(("+data_E+
")/(2*(1+("+data_nu+
"))))";
6285 std::string sigma_d = mu+
"*(Grad_"+varname+
"+Grad_"+varname+
"')";
6286 std::string expr =
"sqrt(3/2)*Norm(Deviator("+sigma_d+
"))";
6287 ga_interpolation_Lagrange_fem(md, expr, mf_vm, VM);
6297 struct linear_incompressibility_brick :
public virtual_brick {
6299 virtual void asm_real_tangent_terms(
const model &md,
size_type ,
6300 const model::varnamelist &vl,
6301 const model::varnamelist &dl,
6302 const model::mimlist &mims,
6303 model::real_matlist &matl,
6304 model::real_veclist &,
6305 model::real_veclist &,
6307 build_version)
const {
6309 GMM_ASSERT1((matl.size() == 1 && dl.size() == 0)
6310 || (matl.size() == 2 && dl.size() == 1),
6311 "Wrong term and/or data number for Linear incompressibility "
6313 GMM_ASSERT1(mims.size() == 1,
"Linear incompressibility brick need one "
6314 "and only one mesh_im");
6315 GMM_ASSERT1(vl.size() == 2,
"Wrong number of variables for linear "
6316 "incompressibility brick");
6318 bool penalized = (dl.size() == 1);
6319 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6320 const mesh_fem &mf_p = md.mesh_fem_of_variable(vl[1]);
6321 const mesh_im &mim = *mims[0];
6322 const model_real_plain_vector *COEFF = 0;
6323 const mesh_fem *mf_data = 0;
6326 COEFF = &(md.real_variable(dl[0]));
6327 mf_data = md.pmesh_fem_of_variable(dl[0]);
6329 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
6330 GMM_ASSERT1(s == 1,
"Bad format for the penalization parameter");
6333 mesh_region rg(region);
6334 mim.linked_mesh().intersect_with_mpi_region(rg);
6336 GMM_TRACE2(
"Stokes term assembly");
6344 gmm::scale(matl[1], scalar_type(-1));
6348 gmm::scale(matl[1], -(*COEFF)[0]);
6355 virtual void real_post_assembly_in_serial(
const model &,
size_type,
6356 const model::varnamelist &,
6357 const model::varnamelist &,
6358 const model::mimlist &,
6359 model::real_matlist &,
6360 model::real_veclist &,
6361 model::real_veclist &,
6363 build_version)
const
6367 linear_incompressibility_brick() {
6368 set_flags(
"Linear incompressibility brick",
6378 const std::string &multname,
size_type region,
6379 const std::string &dataexpr) {
6381 pbrick pbr = std::make_shared<linear_incompressibility_brick>();
6383 tl.push_back(model::term_description(multname, varname,
true));
6384 model::varnamelist vl(1, varname);
6385 vl.push_back(multname);
6386 model::varnamelist dl;
6387 if (dataname.size()) {
6388 dl.push_back(dataexpr);
6389 tl.push_back(model::term_description(multname, multname,
true));
6391 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
6393 std::string test_varname
6394 =
"Test_" + sup_previous_and_dot_to_varname(varname);
6395 std::string test_multname
6396 =
"Test_" + sup_previous_and_dot_to_varname(multname);
6398 if (dataexpr.size())
6399 expr =
"-"+multname+
"*Div_"+test_varname +
"-"+test_multname
6400 +
"*Div_"+varname+
"+(("+dataexpr+
")*"+multname+
")*"+test_multname;
6402 expr =
"-"+multname+
"*Div_"+test_varname +
"-"+test_multname
6405 "Linear incompressibility",
true);
6408 (md, mim, expr, region,
false,
false,
6409 "Linear incompressibility (with nonlinear dependance)");
6422 struct mass_brick :
public virtual_brick {
6424 virtual void asm_real_tangent_terms(
const model &md,
size_type,
6425 const model::varnamelist &vl,
6426 const model::varnamelist &dl,
6427 const model::mimlist &mims,
6428 model::real_matlist &matl,
6429 model::real_veclist &,
6430 model::real_veclist &,
6432 build_version)
const {
6433 GMM_ASSERT1(matl.size() == 1,
6434 "Mass brick has one and only one term");
6435 GMM_ASSERT1(mims.size() == 1,
6436 "Mass brick need one and only one mesh_im");
6437 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
6438 "Wrong number of variables for mass brick");
6440 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6441 const mesh &m = mf_u.linked_mesh();
6442 const mesh_im &mim = *mims[0];
6443 mesh_region rg(region);
6444 m.intersect_with_mpi_region(rg);
6446 const mesh_fem *mf_rho = 0;
6447 const model_real_plain_vector *rho = 0;
6450 mf_rho = md.pmesh_fem_of_variable(dl[0]);
6451 rho = &(md.real_variable(dl[0]));
6453 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6454 GMM_ASSERT1(sl == 1,
"Bad format of mass brick coefficient");
6457 GMM_TRACE2(
"Mass matrix assembly");
6459 if (dl.size() && mf_rho) {
6463 if (dl.size()) gmm::scale(matl[0], (*rho)[0]);
6467 virtual void asm_complex_tangent_terms(
const model &md,
size_type,
6468 const model::varnamelist &vl,
6469 const model::varnamelist &dl,
6470 const model::mimlist &mims,
6471 model::complex_matlist &matl,
6472 model::complex_veclist &,
6473 model::complex_veclist &,
6475 build_version)
const {
6476 GMM_ASSERT1(matl.size() == 1,
6477 "Mass brick has one and only one term");
6478 GMM_ASSERT1(mims.size() == 1,
6479 "Mass brick need one and only one mesh_im");
6480 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
6481 "Wrong number of variables for mass brick");
6483 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6484 const mesh &m = mf_u.linked_mesh();
6485 const mesh_im &mim = *mims[0];
6486 mesh_region rg(region);
6487 m.intersect_with_mpi_region(rg);
6489 const mesh_fem *mf_rho = 0;
6490 const model_complex_plain_vector *rho = 0;
6493 mf_rho = md.pmesh_fem_of_variable(dl[0]);
6494 rho = &(md.complex_variable(dl[0]));
6496 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6497 GMM_ASSERT1(sl == 1,
"Bad format of mass brick coefficient");
6500 GMM_TRACE2(
"Mass matrix assembly");
6502 if (dl.size() && mf_rho) {
6506 if (dl.size()) gmm::scale(matl[0], (*rho)[0]);
6510 virtual std::string declare_volume_assembly_string
6511 (
const model &,
size_type,
const model::varnamelist &,
6512 const model::varnamelist &)
const {
6513 return std::string();
6517 set_flags(
"Mass brick",
true ,
6527 const std::string &dataexpr_rho,
size_type region) {
6529 pbrick pbr = std::make_shared<mass_brick>();
6531 tl.push_back(model::term_description(varname, varname,
true));
6532 model::varnamelist dl;
6533 if (dataexpr_rho.size())
6534 dl.push_back(dataexpr_rho);
6535 return md.
add_brick(pbr, model::varnamelist(1, varname), dl, tl,
6536 model::mimlist(1, &mim), region);
6538 std::string test_varname
6539 =
"Test_" + sup_previous_and_dot_to_varname(varname);
6541 if (dataexpr_rho.size())
6542 expr =
"(("+dataexpr_rho+
")*"+varname+
")."+test_varname;
6544 expr = varname+
"."+test_varname;
6546 "Mass matrix",
true);
6549 "Mass matrix (nonlinear)");
6560 struct lumped_mass_for_first_order_brick :
public virtual_brick {
6562 virtual void asm_real_tangent_terms(
const model &md,
size_type,
6563 const model::varnamelist &vl,
6564 const model::varnamelist &dl,
6565 const model::mimlist &mims,
6566 model::real_matlist &matl,
6567 model::real_veclist &,
6568 model::real_veclist &,
6570 build_version)
const {
6571 GMM_ASSERT1(matl.size() == 1,
6572 "Lumped Mass brick has one and only one term");
6573 GMM_ASSERT1(mims.size() == 1,
6574 "Lumped Mass brick needs one and only one mesh_im");
6575 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
6576 "Wrong number of variables for lumped mass brick");
6578 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6579 const mesh &m = mf_u.linked_mesh();
6580 const mesh_im &mim = *mims[0];
6581 mesh_region rg(region);
6582 m.intersect_with_mpi_region(rg);
6584 const mesh_fem *mf_rho = 0;
6585 const model_real_plain_vector *rho = 0;
6588 mf_rho = md.pmesh_fem_of_variable(dl[0]);
6589 rho = &(md.real_variable(dl[0]));
6591 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6592 GMM_ASSERT1(sl == 1,
"Bad format of mass brick coefficient");
6595 GMM_TRACE2(
"Lumped mass matrix assembly (please check that integration is 1st order.)");
6597 if (dl.size() && mf_rho) {
6601 if (dl.size()) gmm::scale(matl[0], (*rho)[0]);
6606 lumped_mass_for_first_order_brick() {
6607 set_flags(
"Lumped mass brick",
true ,
6617 const std::string &dataexpr_rho,
size_type region) {
6618 pbrick pbr = std::make_shared<lumped_mass_for_first_order_brick>();
6620 tl.push_back(model::term_description(varname, varname,
true));
6621 model::varnamelist dl;
6622 if (dataexpr_rho.size())
6623 dl.push_back(dataexpr_rho);
6624 return md.
add_brick(pbr, model::varnamelist(1, varname), dl, tl,
6625 model::mimlist(1, &mim), region);
6641 struct basic_d_on_dt_brick :
public virtual_brick {
6643 virtual void asm_real_tangent_terms(
const model &md,
size_type ib,
6644 const model::varnamelist &vl,
6645 const model::varnamelist &dl,
6646 const model::mimlist &mims,
6647 model::real_matlist &matl,
6648 model::real_veclist &vecl,
6649 model::real_veclist &,
6651 build_version version)
const {
6652 GMM_ASSERT1(matl.size() == 1,
6653 "Basic d/dt brick has one and only one term");
6654 GMM_ASSERT1(mims.size() == 1,
6655 "Basic d/dt brick need one and only one mesh_im");
6656 GMM_ASSERT1(vl.size() == 1 && dl.size() >= 2 && dl.size() <= 3,
6657 "Wrong number of variables for basic d/dt brick");
6661 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
6662 || (md.is_var_newer_than_brick(dl[1], ib));
6664 recompute_matrix = recompute_matrix ||
6665 md.is_var_newer_than_brick(dl[2], ib);
6668 if (recompute_matrix) {
6669 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6670 const mesh &m = mf_u.linked_mesh();
6671 const mesh_im &mim = *mims[0];
6672 mesh_region rg(region);
6673 m.intersect_with_mpi_region(rg);
6675 const model_real_plain_vector &dt = md.real_variable(dl[1]);
6676 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for time step");
6678 const mesh_fem *mf_rho = 0;
6679 const model_real_plain_vector *rho = 0;
6681 if (dl.size() > 2) {
6682 mf_rho = md.pmesh_fem_of_variable(dl[2]);
6683 rho = &(md.real_variable(dl[2]));
6685 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6686 GMM_ASSERT1(sl == 1,
"Bad format for density");
6689 GMM_TRACE2(
"Mass matrix assembly for d_on_dt brick");
6690 if (dl.size() > 2 && mf_rho) {
6693 gmm::scale(matl[0], scalar_type(1) / dt[0]);
6697 if (dl.size() > 2) gmm::scale(matl[0], (*rho)[0] / dt[0]);
6698 else gmm::scale(matl[0], scalar_type(1) / dt[0]);
6701 gmm::mult(matl[0], md.real_variable(dl[0], 1), vecl[0]);
6704 virtual void asm_complex_tangent_terms(
const model &md,
size_type ib,
6705 const model::varnamelist &vl,
6706 const model::varnamelist &dl,
6707 const model::mimlist &mims,
6708 model::complex_matlist &matl,
6709 model::complex_veclist &vecl,
6710 model::complex_veclist &,
6712 build_version version)
const {
6713 GMM_ASSERT1(matl.size() == 1,
6714 "Basic d/dt brick has one and only one term");
6715 GMM_ASSERT1(mims.size() == 1,
6716 "Basic d/dt brick need one and only one mesh_im");
6717 GMM_ASSERT1(vl.size() == 1 && dl.size() >= 2 && dl.size() <= 3,
6718 "Wrong number of variables for basic d/dt brick");
6722 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
6723 || (md.is_var_newer_than_brick(dl[1], ib));
6725 recompute_matrix = recompute_matrix ||
6726 md.is_var_newer_than_brick(dl[2], ib);
6728 if (recompute_matrix) {
6729 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6730 const mesh &m = mf_u.linked_mesh();
6731 const mesh_im &mim = *mims[0];
6733 mesh_region rg(region);
6734 m.intersect_with_mpi_region(rg);
6736 const model_complex_plain_vector &dt = md.complex_variable(dl[1]);
6737 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for time step");
6739 const mesh_fem *mf_rho = 0;
6740 const model_complex_plain_vector *rho = 0;
6742 if (dl.size() > 2) {
6743 mf_rho = md.pmesh_fem_of_variable(dl[2]);
6744 rho = &(md.complex_variable(dl[2]));
6746 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6747 GMM_ASSERT1(sl == 1,
"Bad format for density");
6750 GMM_TRACE2(
"Mass matrix assembly for d_on_dt brick");
6751 if (dl.size() > 2 && mf_rho) {
6754 gmm::scale(matl[0], scalar_type(1) / dt[0]);
6758 if (dl.size() > 2) gmm::scale(matl[0], (*rho)[0] / dt[0]);
6759 else gmm::scale(matl[0], scalar_type(1) / dt[0]);
6762 gmm::mult(matl[0], md.complex_variable(dl[0], 1), vecl[0]);
6765 virtual std::string declare_volume_assembly_string
6766 (
const model &,
size_type,
const model::varnamelist &,
6767 const model::varnamelist &)
const {
6768 return std::string();
6771 basic_d_on_dt_brick() {
6772 set_flags(
"Basic d/dt brick",
true ,
6782 const std::string &dataname_dt,
const std::string &dataname_rho,
6784 pbrick pbr = std::make_shared<basic_d_on_dt_brick>();
6786 tl.push_back(model::term_description(varname, varname,
true));
6787 model::varnamelist dl(1, varname);
6788 dl.push_back(dataname_dt);
6789 if (dataname_rho.size())
6790 dl.push_back(dataname_rho);
6791 return md.
add_brick(pbr, model::varnamelist(1, varname), dl, tl,
6792 model::mimlist(1, &mim), region);
6803 struct basic_d2_on_dt2_brick :
public virtual_brick {
6805 mutable scalar_type old_alphadt2;
6807 virtual void asm_real_tangent_terms(
const model &md,
size_type ib,
6808 const model::varnamelist &vl,
6809 const model::varnamelist &dl,
6810 const model::mimlist &mims,
6811 model::real_matlist &matl,
6812 model::real_veclist &vecl,
6813 model::real_veclist &,
6815 build_version version)
const {
6816 GMM_ASSERT1(matl.size() == 1,
6817 "Basic d2/dt2 brick has one and only one term");
6818 GMM_ASSERT1(mims.size() == 1,
6819 "Basic d2/dt2 brick need one and only one mesh_im");
6820 GMM_ASSERT1(vl.size() == 1 && dl.size() >= 4 && dl.size() <= 5,
6821 "Wrong number of variables for basic d2/dt2 brick");
6823 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0);
6825 recompute_matrix = recompute_matrix || md.is_var_newer_than_brick(dl[2], ib);
6826 if (dl.size() > 4) recompute_matrix || md.is_var_newer_than_brick(dl[4], ib);
6828 const model_real_plain_vector &dt = md.real_variable(dl[2]);
6829 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for time step");
6830 const model_real_plain_vector &alpha = md.real_variable(dl[3]);
6831 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for parameter alpha");
6832 scalar_type alphadt2 = gmm::sqr(dt[0]) * alpha[0];
6834 if (!recompute_matrix && alphadt2 != old_alphadt2)
6835 gmm::scale(matl[0], old_alphadt2/alphadt2);
6836 old_alphadt2 = alphadt2;
6838 if (recompute_matrix) {
6839 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6840 const mesh &m = mf_u.linked_mesh();
6841 const mesh_im &mim = *mims[0];
6842 mesh_region rg(region);
6843 m.intersect_with_mpi_region(rg);
6845 const mesh_fem *mf_rho = 0;
6846 const model_real_plain_vector *rho = 0;
6848 if (dl.size() > 4) {
6849 mf_rho = md.pmesh_fem_of_variable(dl[4]);
6850 rho = &(md.real_variable(dl[4]));
6852 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6853 GMM_ASSERT1(sl == 1,
"Bad format for density");
6856 GMM_TRACE2(
"Mass matrix assembly for d2_on_dt2 brick");
6857 if (dl.size() > 4 && mf_rho) {
6860 gmm::scale(matl[0], scalar_type(1) / alphadt2);
6864 if (dl.size() > 4) gmm::scale(matl[0], (*rho)[0] / alphadt2);
6865 else gmm::scale(matl[0], scalar_type(1) / alphadt2);
6868 gmm::mult(matl[0], md.real_variable(dl[0], 1), vecl[0]);
6869 gmm::mult_add(matl[0], gmm::scaled(md.real_variable(dl[1], 1), dt[0]),
6873 virtual void asm_complex_tangent_terms(
const model &md,
size_type ib,
6874 const model::varnamelist &vl,
6875 const model::varnamelist &dl,
6876 const model::mimlist &mims,
6877 model::complex_matlist &matl,
6878 model::complex_veclist &vecl,
6879 model::complex_veclist &,
6881 build_version version)
const {
6882 GMM_ASSERT1(matl.size() == 1,
6883 "Basic d2/dt2 brick has one and only one term");
6884 GMM_ASSERT1(mims.size() == 1,
6885 "Basic d2/dt2 brick need one and only one mesh_im");
6886 GMM_ASSERT1(vl.size() == 1 && dl.size() >= 4 && dl.size() <= 5,
6887 "Wrong number of variables for basic d2/dt2 brick");
6889 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0);
6891 recompute_matrix = recompute_matrix || md.is_var_newer_than_brick(dl[2], ib);
6892 if (dl.size() > 4) recompute_matrix || md.is_var_newer_than_brick(dl[4], ib);
6895 const model_complex_plain_vector &dt = md.complex_variable(dl[2]);
6896 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for time step");
6897 const model_complex_plain_vector &
alpha = md.complex_variable(dl[3]);
6898 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for parameter alpha");
6899 scalar_type alphadt2 = gmm::real(gmm::sqr(dt[0]) * alpha[0]);
6901 if (!recompute_matrix && alphadt2 != old_alphadt2)
6902 gmm::scale(matl[0], old_alphadt2/alphadt2);
6903 old_alphadt2 = alphadt2;
6905 if (recompute_matrix) {
6906 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6907 const mesh &m = mf_u.linked_mesh();
6908 const mesh_im &mim = *mims[0];
6909 mesh_region rg(region);
6910 m.intersect_with_mpi_region(rg);
6912 const mesh_fem *mf_rho = 0;
6913 const model_complex_plain_vector *rho = 0;
6915 if (dl.size() > 4) {
6916 mf_rho = md.pmesh_fem_of_variable(dl[4]);
6917 rho = &(md.complex_variable(dl[4]));
6919 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6920 GMM_ASSERT1(sl == 1,
"Bad format for density");
6923 GMM_TRACE2(
"Mass matrix assembly for d2_on_dt2 brick");
6924 if (dl.size() > 4 && mf_rho) {
6927 gmm::scale(matl[0], scalar_type(1) / alphadt2);
6931 if (dl.size() > 4) gmm::scale(matl[0], (*rho)[0] / alphadt2);
6932 else gmm::scale(matl[0], scalar_type(1) / alphadt2);
6935 gmm::mult(matl[0], md.complex_variable(dl[0], 1), vecl[0]);
6936 gmm::mult_add(matl[0], gmm::scaled(md.complex_variable(dl[1], 1), dt[0]),
6940 virtual std::string declare_volume_assembly_string
6941 (
const model &,
size_type,
const model::varnamelist &,
6942 const model::varnamelist &)
const {
6943 return std::string();
6946 basic_d2_on_dt2_brick() {
6947 set_flags(
"Basic d2/dt2 brick",
true ,
6957 const std::string &datanameV,
6958 const std::string &dataname_dt,
6959 const std::string &dataname_alpha,
6960 const std::string &dataname_rho,
6962 pbrick pbr = std::make_shared<basic_d2_on_dt2_brick>();
6964 tl.push_back(model::term_description(varnameU, varnameU,
true));
6965 model::varnamelist dl(1, varnameU);
6966 dl.push_back(datanameV);
6967 dl.push_back(dataname_dt);
6968 dl.push_back(dataname_alpha);
6969 if (dataname_rho.size())
6970 dl.push_back(dataname_rho);
6971 return md.
add_brick(pbr, model::varnamelist(1, varnameU), dl, tl,
6972 model::mimlist(1, &mim), region);
6991 void theta_method_dispatcher::set_dispatch_coeff(
const model &md,
size_type ib)
const {
6993 if (md.is_complex())
6994 theta = gmm::real(md.complex_variable(param_names[0])[0]);
6996 theta = md.real_variable(param_names[0])[0];
6998 md.matrix_coeff_of_brick(ib) = theta;
7000 md.rhs_coeffs_of_brick(ib)[0] = theta;
7002 md.rhs_coeffs_of_brick(ib)[1] = (scalar_type(1) - theta);
7006 void theta_method_dispatcher::next_real_iter
7007 (
const model &md,
size_type ib,
const model::varnamelist &vl,
7008 const model::varnamelist &dl, model::real_matlist &matl,
7009 std::vector<model::real_veclist> &vectl,
7010 std::vector<model::real_veclist> &vectl_sym,
bool first_iter)
const {
7011 next_iter(md, ib, vl, dl, matl, vectl, vectl_sym, first_iter);
7014 void theta_method_dispatcher::next_complex_iter
7015 (
const model &md,
size_type ib,
const model::varnamelist &vl,
7016 const model::varnamelist &dl,
7017 model::complex_matlist &matl,
7018 std::vector<model::complex_veclist> &vectl,
7019 std::vector<model::complex_veclist> &vectl_sym,
7020 bool first_iter)
const {
7021 next_iter(md, ib, vl, dl, matl, vectl, vectl_sym, first_iter);
7024 void theta_method_dispatcher::asm_real_tangent_terms
7025 (
const model &md,
size_type ib, model::real_matlist &,
7026 std::vector<model::real_veclist> &,
7027 std::vector<model::real_veclist> &,
7028 build_version version)
const
7029 { md.brick_call(ib, version, 0); }
7031 void theta_method_dispatcher::asm_complex_tangent_terms
7032 (
const model &md,
size_type ib, model::complex_matlist &,
7033 std::vector<model::complex_veclist> &,
7034 std::vector<model::complex_veclist> &,
7035 build_version version)
const
7036 { md.brick_call(ib, version, 0); }
7038 theta_method_dispatcher::theta_method_dispatcher(
const std::string &THETA)
7039 : virtual_dispatcher(2) {
7040 param_names.push_back(THETA);
7044 (
model &md, dal::bit_vector ibricks,
const std::string &THETA) {
7045 pdispatcher pdispatch = std::make_shared<theta_method_dispatcher>(THETA);
7046 for (dal::bv_visitor i(ibricks); !i.finished(); ++i)
7051 (
model &md,
const std::string &U,
const std::string &V,
7052 const std::string &pdt,
const std::string &ptheta) {
7058 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for time step");
7060 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for parameter theta");
7063 scalar_type(1) - scalar_type(1) / theta[0]),
7066 scalar_type(1) / (theta[0]*dt[0])),
7069 -scalar_type(1) / (theta[0]*dt[0])),
7073 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for time step");
7074 const model_real_plain_vector &theta = md.
real_variable(ptheta);
7075 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for parameter theta");
7078 scalar_type(1) - scalar_type(1) / theta[0]),
7081 scalar_type(1) / (theta[0]*dt[0])),
7084 -scalar_type(1) / (theta[0]*dt[0])),
7097 const std::string &pdt,
const std::string &ptwobeta,
7098 const std::string &pgamma) {
7108 if (twobeta != gamma) {
7110 md.set_dispatch_coeff();
7114 md.
assembly(model::BUILD_RHS_WITH_LIN);
7117 model_complex_plain_vector W(nbdof), RHS(nbdof);
7118 gmm::copy(gmm::sub_vector(md.
complex_rhs(), md.interval_of_variable(U)),
7123 gmm::cg(md.linear_complex_matrix_term(id2dt2b, 0),
7124 W, RHS, gmm::identity_matrix(), iter);
7125 GMM_ASSERT1(iter.converged(),
"Velocity not well computed");
7127 gmm::scaled(W, complex_type(1)/(twobeta*dt)),
7131 if (twobeta != gamma) {
7133 md.set_dispatch_coeff();
7137 GMM_ASSERT1(
false,
"to be done");
7146 if (twobeta != gamma) {
7148 md.set_dispatch_coeff();
7152 md.
assembly(model::BUILD_RHS_WITH_LIN);
7155 model_real_plain_vector W(nbdof), RHS(nbdof);
7156 gmm::copy(gmm::sub_vector(md.
real_rhs(), md.interval_of_variable(U)),
7161 gmm::cg(md.linear_real_matrix_term(id2dt2b, 0),
7162 W, RHS, gmm::identity_matrix(), iter);
7163 GMM_ASSERT1(iter.converged(),
"Velocity not well computed");
7165 gmm::scaled(W, scalar_type(1)/(twobeta*dt)),
7169 if (twobeta != gamma) {
7171 md.set_dispatch_coeff();
7186 class midpoint_dispatcher :
public virtual_dispatcher {
7188 gmm::uint64_type id_num;
7192 typedef model::build_version build_version;
7194 void set_dispatch_coeff(
const model &md,
size_type ib)
const {
7195 md.matrix_coeff_of_brick(ib) = scalar_type(1)/scalar_type(2);
7196 md.rhs_coeffs_of_brick(ib)[0] = scalar_type(1);
7197 md.rhs_coeffs_of_brick(ib)[1] = scalar_type(1)/scalar_type(2);
7200 template <
typename MATLIST,
typename VECTLIST>
7201 inline void next_iter(
const model &md,
size_type ib,
7202 const model::varnamelist &vl,
7203 const model::varnamelist &dl,
7205 VECTLIST &vectl, VECTLIST &vectl_sym,
7206 bool first_iter)
const {
7208 pbrick pbr = md.brick_pointer(ib);
7212 if (!(pbr->is_linear()))
7213 md.add_temporaries(vl, id_num);
7214 md.add_temporaries(dl, id_num);
7219 if (pbr->is_linear()) {
7222 if (first_iter) md.update_brick(ib, model::BUILD_RHS);
7225 md.linear_brick_add_to_rhs(ib, 1, 0);
7230 (
const model &md,
size_type ib,
const model::varnamelist &vl,
7231 const model::varnamelist &dl, model::real_matlist &matl,
7232 std::vector<model::real_veclist> &vectl,
7233 std::vector<model::real_veclist> &vectl_sym,
bool first_iter)
const {
7234 next_iter(md, ib, vl, dl, matl, vectl, vectl_sym, first_iter);
7237 void next_complex_iter
7238 (
const model &md,
size_type ib,
const model::varnamelist &vl,
7239 const model::varnamelist &dl,
7240 model::complex_matlist &matl,
7241 std::vector<model::complex_veclist> &vectl,
7242 std::vector<model::complex_veclist> &vectl_sym,
7243 bool first_iter)
const {
7244 next_iter(md, ib, vl, dl, matl, vectl, vectl_sym, first_iter);
7247 void asm_real_tangent_terms
7248 (
const model &md,
size_type ib, model::real_matlist &,
7249 std::vector<model::real_veclist> &vectl,
7250 std::vector<model::real_veclist> &vectl_sym,
7251 build_version version)
const {
7253 scalar_type half = scalar_type(1)/scalar_type(2);
7254 pbrick pbr = md.brick_pointer(ib);
7257 const model::varnamelist &vl = md.varnamelist_of_brick(ib);
7258 const model::varnamelist &dl = md.datanamelist_of_brick(ib);
7260 if (!(pbr->is_linear())) {
7261 for (
size_type i = 0; i < vl.size(); ++i) {
7262 bool is_uptodate = md.temporary_uptodate(vl[i], id_num, ind);
7263 if (!is_uptodate && ind !=
size_type(-1))
7264 gmm::add(gmm::scaled(md.real_variable(vl[i], 0), half),
7265 gmm::scaled(md.real_variable(vl[i], 1), half),
7266 md.set_real_variable(vl[i], ind));
7267 md.set_default_iter_of_variable(vl[i], ind);
7272 for (
size_type i = 0; i < dl.size(); ++i) {
7273 bool is_uptodate = md.temporary_uptodate(dl[i], id_num, ind);
7274 if (!is_uptodate && ind !=
size_type(-1)) {
7275 gmm::add(gmm::scaled(md.real_variable(dl[i], 0), half),
7276 gmm::scaled(md.real_variable(dl[i], 1), half),
7277 md.set_real_variable(dl[i], ind));
7279 md.set_default_iter_of_variable(dl[i], ind);
7283 md.brick_call(ib, version, 0);
7284 if (pbr->is_linear()) {
7288 md.linear_brick_add_to_rhs(ib, 1, 1);
7291 md.reset_default_iter_of_variables(dl);
7292 if (!(pbr->is_linear()))
7293 md.reset_default_iter_of_variables(vl);
7296 virtual void asm_complex_tangent_terms
7297 (
const model &md,
size_type ib, model::complex_matlist &,
7298 std::vector<model::complex_veclist> &vectl,
7299 std::vector<model::complex_veclist> &vectl_sym,
7300 build_version version)
const {
7302 scalar_type half = scalar_type(1)/scalar_type(2);
7303 pbrick pbr = md.brick_pointer(ib);
7306 const model::varnamelist &vl = md.varnamelist_of_brick(ib);
7307 const model::varnamelist &dl = md.datanamelist_of_brick(ib);
7309 if (!(pbr->is_linear())) {
7310 for (
size_type i = 0; i < vl.size(); ++i) {
7311 bool is_uptodate = md.temporary_uptodate(vl[i], id_num, ind);
7312 if (!is_uptodate && ind !=
size_type(-1))
7313 gmm::add(gmm::scaled(md.complex_variable(vl[i], 0), half),
7314 gmm::scaled(md.complex_variable(vl[i], 1), half),
7315 md.set_complex_variable(vl[i], ind));
7316 md.set_default_iter_of_variable(vl[i], ind);
7321 for (
size_type i = 0; i < dl.size(); ++i) {
7322 bool is_uptodate = md.temporary_uptodate(dl[i], id_num, ind);
7323 if (!is_uptodate && ind !=
size_type(-1)) {
7324 gmm::add(gmm::scaled(md.complex_variable(dl[i], 0), half),
7325 gmm::scaled(md.complex_variable(dl[i], 1), half),
7326 md.set_complex_variable(dl[i], ind));
7328 md.set_default_iter_of_variable(dl[i], ind);
7332 md.brick_call(ib, version, 0);
7333 if (pbr->is_linear()) {
7337 md.linear_brick_add_to_rhs(ib, 1, 1);
7340 md.reset_default_iter_of_variables(dl);
7341 if (!(pbr->is_linear()))
7342 md.reset_default_iter_of_variables(vl);
7345 midpoint_dispatcher() : virtual_dispatcher(2)
7346 { id_num = act_counter(); }
7351 pdispatcher pdispatch = std::make_shared<midpoint_dispatcher>();
7352 for (dal::bv_visitor i(ibricks); !i.finished(); ++i)