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 = mf ? 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.mf && !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.mf,
"The primal variable of the "
395 "multiplier should be 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.mf && !(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);
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);
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;
593 gmm::range_basis(MM, columns);
594 if (columns.size() == 0)
595 GMM_WARNING1(
"Empty basis found for multiplier " << multname);
597 if (mults.size() > 1) {
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) {
622 gmm::range_basis(MGLOB, glob_columns, 1E-12, gmm::col_major(),
true);
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.mf ?
"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].mf == 0,
755 "Cannot explicitly resize a fem variable or data");
756 GMM_ASSERT1(variables[name].imd == 0,
757 "Cannot explicitly resize an im variable or 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.mf && (v.second.filter & VDESCRFILTER_INFSUP) &&
973 mim == v.second.filter_mim) found =
true;
975 if (!found) sup_dependency(*mim);
978 is_linear_ = is_symmetric_ = is_coercive_ =
true;
979 for (dal::bv_visitor ibb(valid_bricks); !ibb.finished(); ++ibb) {
980 is_linear_ = is_linear_ && bricks[ibb].pbr->is_linear();
981 is_symmetric_ = is_symmetric_ && bricks[ibb].pbr->is_symmetric();
982 is_coercive_ = is_coercive_ && bricks[ibb].pbr->is_coercive();
984 bricks[ib] = brick_description();
988 for (dal::bv_visitor ibb(valid_bricks); !ibb.finished(); ++ibb) {
989 for (
const auto &vname : bricks[ibb].vlist)
990 GMM_ASSERT1(varname.compare(vname),
991 "Cannot delete a variable which is still used by a brick");
992 for (
const auto &dname : bricks[ibb].dlist)
993 GMM_ASSERT1(varname.compare(dname),
994 "Cannot delete a data which is still used by a brick");
997 VAR_SET::const_iterator it = find_variable(varname);
1000 const mesh_fem *mf = it->second.mf;
1002 for(VAR_SET::iterator it2 = variables.begin();
1003 it2 != variables.end(); ++it2) {
1004 if (it != it2 && it2->second.mf && mf == it2->second.mf)
1007 if (!found) sup_dependency(*mf);
1009 if (it->second.filter & VDESCRFILTER_INFSUP) {
1010 const mesh_im *mim = it->second.filter_mim;
1012 for (dal::bv_visitor ibb(valid_bricks); !ibb.finished(); ++ibb) {
1013 for (
size_type j = 0; j < bricks[ibb].mims.size(); ++j)
1014 if (bricks[ibb].mims[j] == mim) found =
true;
1016 for (VAR_SET::iterator it2 = variables.begin();
1017 it2 != variables.end(); ++it2) {
1018 if (it != it2 && it2->second.mf &&
1019 (it2->second.filter & VDESCRFILTER_INFSUP) &&
1020 mim == it2->second.filter_mim) found =
true;
1022 if (!found) sup_dependency(*mim);
1026 if (it->second.imd != 0) sup_dependency(*(it->second.imd));
1028 variables.erase(varname);
1029 act_size_to_be_done =
true;
1033 const varnamelist &datanames,
1034 const termlist &terms,
1035 const mimlist &mims,
size_type region) {
1036 size_type ib = valid_bricks.first_false();
1038 for (
size_type i = 0; i < terms.size(); ++i)
1039 if (terms[i].is_global && terms[i].is_matrix_term && pbr->is_linear())
1040 GMM_ASSERT1(
false,
"Global linear matrix terms are not allowed");
1042 if (ib == bricks.size())
1043 bricks.push_back(brick_description(pbr, varnames, datanames, terms,
1046 bricks[ib] = brick_description(pbr, varnames, datanames, terms,
1048 active_bricks.add(ib);
1049 valid_bricks.add(ib);
1056 "Impossible to add a complex brick to a real model");
1058 bricks[ib].cmatlist.resize(terms.size());
1059 bricks[ib].cveclist[0].resize(terms.size());
1060 bricks[ib].cveclist_sym[0].resize(terms.size());
1062 bricks[ib].rmatlist.resize(terms.size());
1063 bricks[ib].rveclist[0].resize(terms.size());
1064 bricks[ib].rveclist_sym[0].resize(terms.size());
1066 is_linear_ = is_linear_ && pbr->is_linear();
1067 is_symmetric_ = is_symmetric_ && pbr->is_symmetric();
1068 is_coercive_ = is_coercive_ && pbr->is_coercive();
1070 for (
const auto &vname : varnames)
1071 GMM_ASSERT1(variables.count(vname),
1072 "Undefined model variable " << vname);
1074 for (
const auto &dname : datanames)
1075 GMM_ASSERT1(variables.count(dname),
1076 "Undefined model data or variable " << dname);
1082 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
1084 bricks[ib].mims.push_back(&mim);
1085 add_dependency(mim);
1089 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
1091 bricks[ib].tlist = terms;
1092 if (
is_complex() && bricks[ib].pbr->is_complex()) {
1093 bricks.back().cmatlist.resize(terms.size());
1094 bricks.back().cveclist[0].resize(terms.size());
1095 bricks.back().cveclist_sym[0].resize(terms.size());
1097 bricks.back().rmatlist.resize(terms.size());
1098 bricks.back().rveclist[0].resize(terms.size());
1099 bricks.back().rveclist_sym[0].resize(terms.size());
1104 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
1106 bricks[ib].vlist = vl;
1107 for (
const auto &v : vl)
1108 GMM_ASSERT1(variables.count(v),
"Undefined model variable " << v);
1112 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
1114 bricks[ib].dlist = dl;
1115 for (
const auto &v : dl)
1116 GMM_ASSERT1(variables.count(v),
"Undefined model variable " << v);
1120 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
1122 bricks[ib].mims = ml;
1123 for (
const auto &mim : ml) add_dependency(*mim);
1127 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
1129 bricks[ib].is_update_brick = flag;
1132 void model::set_time(scalar_type t,
bool to_init) {
1133 static const std::string varname(
"t");
1134 VAR_SET::iterator it = variables.find(varname);
1135 if (it == variables.end()) {
1138 GMM_ASSERT1(it->second.size() == 1,
"Time data should be of size 1");
1140 if (it == variables.end() || to_init) {
1148 scalar_type model::get_time() {
1149 static const std::string varname(
"t");
1150 set_time(scalar_type(0),
false);
1157 void model::call_init_affine_dependent_variables(
int version) {
1158 for (VAR_SET::iterator it = variables.begin();
1159 it != variables.end(); ++it) {
1160 var_description &vdescr = it->second;
1161 if (vdescr.is_variable && vdescr.ptsc) {
1163 vdescr.ptsc->init_affine_dependent_variables_precomputation(*
this);
1165 vdescr.ptsc->init_affine_dependent_variables(*
this);
1170 void model::shift_variables_for_time_integration() {
1171 for (VAR_SET::iterator it = variables.begin();
1172 it != variables.end(); ++it)
1173 if (it->second.is_variable && it->second.ptsc)
1174 it->second.ptsc->shift_variables(*
this);
1177 void model::add_time_integration_scheme(
const std::string &varname,
1178 ptime_scheme ptsc) {
1179 VAR_SET::iterator it = variables.find(varname);
1180 GMM_ASSERT1(it != variables.end(),
"Undefined variable " << varname);
1181 GMM_ASSERT1(it->second.is_variable && !(it->second.is_affine_dependent),
1182 "Cannot apply an integration scheme to " << varname);
1183 it->second.ptsc = ptsc;
1190 time_integration = 1;
1193 void model::copy_init_time_derivative() {
1195 for (VAR_SET::iterator it = variables.begin();
1196 it != variables.end(); ++it)
1197 if (it->second.is_variable && it->second.ptsc) {
1199 std::string name_v, name_previous_v;
1200 it->second.ptsc->time_derivative_to_be_initialized(name_v,
1203 if (name_v.size()) {
1221 class APIDECL first_order_theta_method_scheme
1222 :
public virtual_time_scheme {
1224 std::string U, U0, V, V0;
1229 virtual void init_affine_dependent_variables(model &md)
const {
1230 scalar_type dt = md.get_time_step();
1231 scalar_type a = scalar_type(1)/(theta*dt);
1232 scalar_type b = (scalar_type(1)-theta)/theta;
1233 md.set_factor_of_variable(V, a);
1234 if (md.is_complex()) {
1235 gmm::add(gmm::scaled(md.complex_variable(U0), -complex_type(a)),
1236 gmm::scaled(md.complex_variable(V0), -complex_type(b)),
1237 md.set_complex_constant_part(V));
1240 gmm::add(gmm::scaled(md.real_variable(U0), -a),
1241 gmm::scaled(md.real_variable(V0), -b),
1242 md.set_real_constant_part(V));
1247 virtual void init_affine_dependent_variables_precomputation(model &md)
1249 scalar_type dt = md.get_time_step();
1250 md.set_factor_of_variable(V, scalar_type(1)/dt);
1251 if (md.is_complex()) {
1252 gmm::copy(gmm::scaled(md.complex_variable(U0), -complex_type(1)/dt),
1253 md.set_complex_constant_part(V));
1256 gmm::copy(gmm::scaled(md.real_variable(U0), -scalar_type(1)/dt),
1257 md.set_real_constant_part(V));
1261 virtual void time_derivative_to_be_initialized
1262 (std::string &name_v, std::string &name_previous_v)
const
1263 {
if (theta != scalar_type(1)) { name_v = V; name_previous_v = V0; } }
1265 virtual void shift_variables(model &md)
const {
1266 if (md.is_complex()) {
1267 gmm::copy(md.complex_variable(U), md.set_complex_variable(U0));
1268 gmm::copy(md.complex_variable(V), md.set_complex_variable(V0));
1270 gmm::copy(md.real_variable(U), md.set_real_variable(U0));
1271 gmm::copy(md.real_variable(V), md.set_real_variable(V0));
1276 first_order_theta_method_scheme(model &md, std::string varname,
1279 U0 =
"Previous_" + U;
1281 V0 =
"Previous_Dot_" + U;
1283 GMM_ASSERT1(theta > scalar_type(0) && theta <= scalar_type(1),
1284 "Invalid value of theta parameter for the theta-method");
1286 if (!(md.variable_exists(V)))
1287 md.add_affine_dependent_variable(V, U);
1288 const mesh_fem *mf = md.pmesh_fem_of_variable(U);
1289 size_type s = md.is_complex() ? gmm::vect_size(md.complex_variable(U))
1290 : gmm::vect_size(md.real_variable(U));
1293 if (!(md.variable_exists(U0))) md.add_fem_data(U0, *mf);
1294 if (!(md.variable_exists(V0))) md.add_fem_data(V0, *mf);
1296 if (!(md.variable_exists(U0))) md.add_fixed_size_data(U0, s);
1297 if (!(md.variable_exists(V0))) md.add_fixed_size_data(V0, s);
1304 void add_theta_method_for_first_order(model &md,
const std::string &varname,
1305 scalar_type theta) {
1307 = std::make_shared<first_order_theta_method_scheme>(md, varname,theta);
1308 md.add_time_integration_scheme(varname, ptsc);
1317 class APIDECL second_order_theta_method_scheme
1318 :
public virtual_time_scheme {
1320 std::string U, U0, V, V0, A, A0;
1326 virtual void init_affine_dependent_variables(model &md)
const {
1327 scalar_type dt = md.get_time_step();
1328 md.set_factor_of_variable(V, scalar_type(1)/(theta*dt));
1329 md.set_factor_of_variable(A, scalar_type(1)/(theta*theta*dt*dt));
1330 if (md.is_complex()) {
1331 gmm::add(gmm::scaled(md.complex_variable(U0),
1332 -complex_type(1)/(theta*dt)),
1333 gmm::scaled(md.complex_variable(V0),
1334 -(complex_type(1)-complex_type(theta))/theta),
1335 md.set_complex_constant_part(V));
1336 gmm::add(gmm::scaled(md.complex_variable(U0),
1337 -complex_type(1)/(theta*theta*dt*dt)),
1338 gmm::scaled(md.complex_variable(A0),
1339 -(complex_type(1)-complex_type(theta))/theta),
1340 md.set_complex_constant_part(A));
1341 gmm::add(gmm::scaled(md.complex_variable(V0),
1342 -complex_type(1)/(theta*theta*dt)),
1343 md.set_complex_constant_part(A));
1347 gmm::add(gmm::scaled(md.real_variable(U0),
1348 -scalar_type(1)/(theta*dt)),
1349 gmm::scaled(md.real_variable(V0),
1350 -(scalar_type(1)-theta)/theta),
1351 md.set_real_constant_part(V));
1352 gmm::add(gmm::scaled(md.real_variable(U0),
1353 -scalar_type(1)/(theta*theta*dt*dt)),
1354 gmm::scaled(md.real_variable(A0),
1355 -(scalar_type(1)-theta)/theta),
1356 md.set_real_constant_part(A));
1357 gmm::add(gmm::scaled(md.real_variable(V0),
1358 -scalar_type(1)/(theta*theta*dt)),
1359 md.set_real_constant_part(A));
1366 virtual void init_affine_dependent_variables_precomputation(model &md)
1368 scalar_type dt = md.get_time_step();
1369 md.set_factor_of_variable(V, scalar_type(1)/dt);
1370 md.set_factor_of_variable(A, scalar_type(1)/(dt*dt));
1371 if (md.is_complex()) {
1372 gmm::copy(gmm::scaled(md.complex_variable(U0),
1373 -complex_type(1)/dt),
1374 md.set_complex_constant_part(V));
1375 gmm::add(gmm::scaled(md.complex_variable(U0),
1376 -complex_type(1)/(dt*dt)),
1377 gmm::scaled(md.complex_variable(V0),
1378 -complex_type(1)/dt),
1379 md.set_complex_constant_part(A));
1381 gmm::copy(gmm::scaled(md.real_variable(U0),
1382 -scalar_type(1)/dt),
1383 md.set_real_constant_part(V));
1384 gmm::add(gmm::scaled(md.real_variable(U0),
1385 -scalar_type(1)/(dt*dt)),
1386 gmm::scaled(md.real_variable(V0),
1387 -scalar_type(1)/dt),
1388 md.set_real_constant_part(A));
1392 virtual void time_derivative_to_be_initialized
1393 (std::string &name_v, std::string &name_previous_v)
const
1394 {
if (theta != scalar_type(1)) { name_v = A; name_previous_v = A0; } }
1396 virtual void shift_variables(model &md)
const {
1397 if (md.is_complex()) {
1398 gmm::copy(md.complex_variable(U), md.set_complex_variable(U0));
1399 gmm::copy(md.complex_variable(V), md.set_complex_variable(V0));
1400 gmm::copy(md.complex_variable(A), md.set_complex_variable(A0));
1402 gmm::copy(md.real_variable(U), md.set_real_variable(U0));
1403 gmm::copy(md.real_variable(V), md.set_real_variable(V0));
1404 gmm::copy(md.real_variable(A), md.set_real_variable(A0));
1409 second_order_theta_method_scheme(model &md, std::string varname,
1412 U0 =
"Previous_" + U;
1414 V0 =
"Previous_Dot_" + U;
1416 A0 =
"Previous_Dot2_" + U;
1418 GMM_ASSERT1(theta > scalar_type(0) && theta <= scalar_type(1),
1419 "Invalid value of theta parameter for the theta-method");
1421 if (!(md.variable_exists(V)))
1422 md.add_affine_dependent_variable(V, U);
1423 if (!(md.variable_exists(A)))
1424 md.add_affine_dependent_variable(A, U);
1426 const mesh_fem *mf = md.pmesh_fem_of_variable(U);
1427 size_type s = md.is_complex() ? gmm::vect_size(md.complex_variable(U))
1428 : gmm::vect_size(md.real_variable(U));
1431 if (!(md.variable_exists(U0))) md.add_fem_data(U0, *mf);
1432 if (!(md.variable_exists(V0))) md.add_fem_data(V0, *mf);
1433 if (!(md.variable_exists(A0))) md.add_fem_data(A0, *mf);
1435 if (!(md.variable_exists(U0))) md.add_fixed_size_data(U0, s);
1436 if (!(md.variable_exists(V0))) md.add_fixed_size_data(V0, s);
1437 if (!(md.variable_exists(A0))) md.add_fixed_size_data(A0, s);
1444 void add_theta_method_for_second_order(model &md,
const std::string &varname,
1445 scalar_type theta) {
1446 ptime_scheme ptsc = std::make_shared<second_order_theta_method_scheme>
1448 md.add_time_integration_scheme(varname, ptsc);
1458 class APIDECL Newmark_scheme
1459 :
public virtual_time_scheme {
1461 std::string U, U0, V, V0, A, A0;
1462 scalar_type beta, gamma;
1467 virtual void init_affine_dependent_variables(model &md)
const {
1468 scalar_type dt = md.get_time_step();
1469 scalar_type a0 = scalar_type(1)/(beta*dt*dt), a1 = dt*a0;
1470 scalar_type a2 = (scalar_type(1) - scalar_type(2)*beta)
1471 / (scalar_type(2)*beta);
1472 scalar_type b0 = gamma/(beta*dt), b1 = (beta-gamma)/beta;
1473 scalar_type b2 = dt*(scalar_type(1)-gamma/(scalar_type(2)*beta));
1475 md.set_factor_of_variable(V, b0);
1476 md.set_factor_of_variable(A, a0);
1477 if (md.is_complex()) {
1478 gmm::add(gmm::scaled(md.complex_variable(U0), -complex_type(b0)),
1479 gmm::scaled(md.complex_variable(V0), complex_type(b1)),
1480 md.set_complex_constant_part(V));
1481 gmm::add(gmm::scaled(md.complex_variable(A0), complex_type(b2)),
1482 md.set_complex_constant_part(V));
1483 gmm::add(gmm::scaled(md.complex_variable(U0), -complex_type(a0)),
1484 gmm::scaled(md.complex_variable(V0), -complex_type(a1)),
1485 md.set_complex_constant_part(A));
1486 gmm::add(gmm::scaled(md.complex_variable(A0), -complex_type(a2)),
1487 md.set_complex_constant_part(A));
1489 gmm::add(gmm::scaled(md.real_variable(U0), -b0),
1490 gmm::scaled(md.real_variable(V0), b1),
1491 md.set_real_constant_part(V));
1492 gmm::add(gmm::scaled(md.real_variable(A0), b2),
1493 md.set_real_constant_part(V));
1494 gmm::add(gmm::scaled(md.real_variable(U0), -a0),
1495 gmm::scaled(md.real_variable(V0), -a1),
1496 md.set_real_constant_part(A));
1497 gmm::add(gmm::scaled(md.real_variable(A0), -a2),
1498 md.set_real_constant_part(A));
1505 virtual void init_affine_dependent_variables_precomputation(model &md)
1507 scalar_type dt = md.get_time_step();
1508 md.set_factor_of_variable(V, scalar_type(1)/dt);
1509 md.set_factor_of_variable(A, scalar_type(1)/(dt*dt));
1510 if (md.is_complex()) {
1511 gmm::copy(gmm::scaled(md.complex_variable(U0),
1512 -complex_type(1)/dt),
1513 md.set_complex_constant_part(V));
1514 gmm::add(gmm::scaled(md.complex_variable(U0),
1515 -complex_type(1)/(dt*dt)),
1516 gmm::scaled(md.complex_variable(V0),
1517 -complex_type(1)/dt),
1518 md.set_complex_constant_part(A));
1520 gmm::copy(gmm::scaled(md.real_variable(U0),
1521 -scalar_type(1)/dt),
1522 md.set_real_constant_part(V));
1523 gmm::add(gmm::scaled(md.real_variable(U0),
1524 -scalar_type(1)/(dt*dt)),
1525 gmm::scaled(md.real_variable(V0),
1526 -scalar_type(1)/dt),
1527 md.set_real_constant_part(A));
1531 virtual void time_derivative_to_be_initialized
1532 (std::string &name_v, std::string &name_previous_v)
const {
1533 if (beta != scalar_type(0.5) || gamma != scalar_type(1))
1534 { name_v = A; name_previous_v = A0; }
1537 virtual void shift_variables(model &md)
const {
1538 if (md.is_complex()) {
1539 gmm::copy(md.complex_variable(U), md.set_complex_variable(U0));
1540 gmm::copy(md.complex_variable(V), md.set_complex_variable(V0));
1541 gmm::copy(md.complex_variable(A), md.set_complex_variable(A0));
1543 gmm::copy(md.real_variable(U), md.set_real_variable(U0));
1544 gmm::copy(md.real_variable(V), md.set_real_variable(V0));
1545 gmm::copy(md.real_variable(A), md.set_real_variable(A0));
1550 Newmark_scheme(model &md, std::string varname,
1551 scalar_type be, scalar_type ga) {
1553 U0 =
"Previous_" + U;
1555 V0 =
"Previous_Dot_" + U;
1557 A0 =
"Previous_Dot2_" + U;
1558 beta = be; gamma = ga;
1559 GMM_ASSERT1(beta > scalar_type(0) && beta <= scalar_type(1)
1560 && gamma >= scalar_type(0.5) && gamma <= scalar_type(1),
1561 "Invalid parameter values for the Newmark scheme");
1563 if (!(md.variable_exists(V)))
1564 md.add_affine_dependent_variable(V, U);
1565 if (!(md.variable_exists(A)))
1566 md.add_affine_dependent_variable(A, U);
1568 const mesh_fem *mf = md.pmesh_fem_of_variable(U);
1569 size_type s = md.is_complex() ? gmm::vect_size(md.complex_variable(U))
1570 : gmm::vect_size(md.real_variable(U));
1573 if (!(md.variable_exists(U0))) md.add_fem_data(U0, *mf);
1574 if (!(md.variable_exists(V0))) md.add_fem_data(V0, *mf);
1575 if (!(md.variable_exists(A0))) md.add_fem_data(A0, *mf);
1577 if (!(md.variable_exists(U0))) md.add_fixed_size_data(U0, s);
1578 if (!(md.variable_exists(V0))) md.add_fixed_size_data(V0, s);
1579 if (!(md.variable_exists(A0))) md.add_fixed_size_data(A0, s);
1586 void add_Newmark_scheme(model &md,
const std::string &varname,
1587 scalar_type beta, scalar_type gamma) {
1588 ptime_scheme ptsc = std::make_shared<Newmark_scheme>
1589 (md, varname, beta, gamma);
1590 md.add_time_integration_scheme(varname, ptsc);
1599 class APIDECL Houbolt_scheme
1600 :
public virtual_time_scheme {
1602 std::string U, U01, U02, U03, V, A;
1607 virtual void init_affine_dependent_variables(model &md)
const {
1608 scalar_type dt = md.get_time_step();
1609 scalar_type a0 = scalar_type(2)/(dt*dt);
1610 scalar_type a1 = scalar_type(5)/(dt*dt);
1611 scalar_type a2 = scalar_type(4)/(dt*dt);
1612 scalar_type a3 = scalar_type(1)/(dt*dt);
1613 scalar_type b0 = scalar_type(11)/(scalar_type(6)*dt);
1614 scalar_type b1 = scalar_type(18)/(scalar_type(6)*dt);
1615 scalar_type b2 = scalar_type(9)/(scalar_type(6)*dt);
1616 scalar_type b3 = scalar_type(2)/(scalar_type(6)*dt);
1618 md.set_factor_of_variable(V, b0);
1619 md.set_factor_of_variable(A, a0);
1620 if (md.is_complex()) {
1621 gmm::add(gmm::scaled(md.complex_variable(U01), -complex_type(b1)),
1622 gmm::scaled(md.complex_variable(U02), complex_type(b2)),
1623 md.set_complex_constant_part(V));
1624 gmm::add(gmm::scaled(md.complex_variable(U03), -complex_type(b3)),
1625 md.set_complex_constant_part(V));
1626 gmm::add(gmm::scaled(md.complex_variable(U01), -complex_type(a1)),
1627 gmm::scaled(md.complex_variable(U02), complex_type(a2)),
1628 md.set_complex_constant_part(A));
1629 gmm::add(gmm::scaled(md.complex_variable(U03), -complex_type(a3)),
1630 md.set_complex_constant_part(A));
1632 gmm::add(gmm::scaled(md.real_variable(U01), -b1),
1633 gmm::scaled(md.real_variable(U02), b2),
1634 md.set_real_constant_part(V));
1635 gmm::add(gmm::scaled(md.real_variable(U03), -b3),
1636 md.set_real_constant_part(V));
1637 gmm::add(gmm::scaled(md.real_variable(U01), -a1),
1638 gmm::scaled(md.real_variable(U02), a2),
1639 md.set_real_constant_part(A));
1640 gmm::add(gmm::scaled(md.real_variable(U03), -a3),
1641 md.set_real_constant_part(A));
1645 virtual void init_affine_dependent_variables_precomputation(model &md)
1650 virtual void time_derivative_to_be_initialized
1651 (std::string &name_v, std::string &name_previous_v)
const {
1653 (void) name_previous_v;
1656 virtual void shift_variables(model &md)
const {
1657 if (md.is_complex()) {
1658 gmm::copy(md.complex_variable(U02), md.set_complex_variable(U03));
1659 gmm::copy(md.complex_variable(U01), md.set_complex_variable(U02));
1660 gmm::copy(md.complex_variable(U), md.set_complex_variable(U01));
1662 gmm::copy(md.real_variable(U02), md.set_real_variable(U03));
1663 gmm::copy(md.real_variable(U01), md.set_real_variable(U02));
1664 gmm::copy(md.real_variable(U), md.set_real_variable(U01));
1669 Houbolt_scheme(model &md, std::string varname) {
1671 U01 =
"Previous_" + U;
1672 U02 =
"Previous2_" + U;
1673 U03 =
"Previous3_" + U;
1677 if (!(md.variable_exists(V)))
1678 md.add_affine_dependent_variable(V, U);
1679 if (!(md.variable_exists(A)))
1680 md.add_affine_dependent_variable(A, U);
1682 const mesh_fem *mf = md.pmesh_fem_of_variable(U);
1683 size_type s = md.is_complex() ? gmm::vect_size(md.complex_variable(U))
1684 : gmm::vect_size(md.real_variable(U));
1687 if (!(md.variable_exists(U01))) md.add_fem_data(U01, *mf);
1688 if (!(md.variable_exists(U02))) md.add_fem_data(U02, *mf);
1689 if (!(md.variable_exists(U03))) md.add_fem_data(U03, *mf);
1691 if (!(md.variable_exists(U01))) md.add_fixed_size_data(U01, s);
1692 if (!(md.variable_exists(U02))) md.add_fixed_size_data(U02, s);
1693 if (!(md.variable_exists(U03))) md.add_fixed_size_data(U03, s);
1700 void add_Houbolt_scheme(model &md,
const std::string &varname) {
1701 ptime_scheme ptsc = std::make_shared<Houbolt_scheme>
1703 md.add_time_integration_scheme(varname, ptsc);
1715 GMM_ASSERT1(valid_bricks[ibrick],
"Inexistent brick");
1716 pbrick pbr = bricks[ibrick].pbr;
1718 bricks[ibrick].pdispatch = pdispatch;
1721 = std::max(
size_type(1), pdispatch->nbrhs());
1723 gmm::resize(bricks[ibrick].coeffs, nbrhs);
1726 bricks[ibrick].cveclist.resize(nbrhs);
1727 bricks[ibrick].cveclist_sym.resize(nbrhs);
1729 bricks[ibrick].cveclist[k] = bricks[ibrick].cveclist[0];
1730 bricks[ibrick].cveclist_sym[k] = bricks[ibrick].cveclist_sym[0];
1733 bricks[ibrick].rveclist.resize(nbrhs);
1734 bricks[ibrick].rveclist_sym.resize(nbrhs);
1736 bricks[ibrick].rveclist[k] = bricks[ibrick].rveclist[0];
1737 bricks[ibrick].rveclist_sym[k] = bricks[ibrick].rveclist_sym[0];
1744 GMM_ASSERT1(valid_bricks[ind_brick],
"Inexistent brick");
1745 GMM_ASSERT1(ind_var < bricks[ind_brick].vlist.size(),
1746 "Inexistent brick variable");
1747 return bricks[ind_brick].vlist[ind_var];
1752 GMM_ASSERT1(valid_bricks[ind_brick],
"Inexistent brick");
1753 GMM_ASSERT1(ind_data < bricks[ind_brick].dlist.size(),
1754 "Inexistent brick data");
1755 return bricks[ind_brick].dlist[ind_data];
1759 if (valid_bricks.card() == 0)
1760 ost <<
"Model with no bricks" << endl;
1762 ost <<
"List of model bricks:" << endl;
1763 for (dal::bv_visitor i(valid_bricks); !i.finished(); ++i) {
1764 ost <<
"Brick " << std::setw(3) << std::right << i + base_id
1765 <<
" " << std::setw(20) << std::right
1766 << bricks[i].pbr->brick_name();
1767 if (!(active_bricks[i])) ost <<
" (deactivated)";
1768 if (bricks[i].pdispatch) ost <<
" (dispatched)";
1769 ost << endl <<
" concerned variables: " << bricks[i].vlist[0];
1770 for (
size_type j = 1; j < bricks[i].vlist.size(); ++j)
1771 ost <<
", " << bricks[i].vlist[j];
1773 ost <<
" brick with " << bricks[i].tlist.size() <<
" term";
1774 if (bricks[i].tlist.size() > 1) ost <<
"s";
1783 void model::brick_init(
size_type ib, build_version version,
1785 const brick_description &brick = bricks[ib];
1786 bool cplx =
is_complex() && brick.pbr->is_complex();
1789 for (
size_type j = 0; j < brick.tlist.size(); ++j) {
1790 const term_description &term = brick.tlist[j];
1791 bool isg = term.is_global;
1793 gmm::vect_size(crhs) : gmm::vect_size(rrhs);
1794 size_type nbd1 = isg ? nbgdof : variables[term.var1].size();
1795 size_type nbd2 = isg ? nbgdof : (term.is_matrix_term ?
1796 variables[term.var2].size() : 0);
1797 if (term.is_matrix_term &&
1798 (brick.pbr->is_linear() || (version | BUILD_MATRIX))) {
1799 if (version | BUILD_ON_DATA_CHANGE) {
1801 gmm::resize(brick.cmatlist[j], nbd1, nbd2);
1803 gmm::resize(brick.rmatlist[j], nbd1, nbd2);
1806 brick.cmatlist[j] = model_complex_sparse_matrix(nbd1, nbd2);
1808 brick.rmatlist[j] = model_real_sparse_matrix(nbd1, nbd2);
1811 if (brick.pbr->is_linear() || (version | BUILD_RHS)) {
1812 for (
size_type k = 0; k < brick.nbrhs; ++k) {
1814 if (k == rhs_ind)
gmm::clear(brick.cveclist[k][j]);
1816 if (term.is_symmetric && term.var1.compare(term.var2)) {
1817 if (k == rhs_ind)
gmm::clear(brick.cveclist_sym[k][j]);
1821 if (k == rhs_ind)
gmm::clear(brick.rveclist[k][j]);
1823 if (term.is_symmetric && term.var1.compare(term.var2)) {
1824 if (k == rhs_ind)
gmm::clear(brick.rveclist_sym[k][j]);
1833 void model::post_to_variables_step(){}
1835 void model::brick_call(
size_type ib, build_version version,
1838 const brick_description &brick = bricks[ib];
1839 bool cplx =
is_complex() && brick.pbr->is_complex();
1841 brick_init(ib, version, rhs_ind);
1845 brick.pbr->complex_pre_assembly_in_serial(*
this, ib, brick.vlist,
1846 brick.dlist, brick.mims,
1848 brick.cveclist[rhs_ind],
1849 brick.cveclist_sym[rhs_ind],
1850 brick.region, version);
1855 accumulated_distro<complex_matlist> cmatlist(brick.cmatlist);
1856 accumulated_distro<complex_veclist> cveclist(brick.cveclist[rhs_ind]);
1857 accumulated_distro<complex_veclist> cveclist_sym(brick.cveclist_sym[rhs_ind]);
1861 brick.pbr->asm_complex_tangent_terms(*
this, ib, brick.vlist,
1862 brick.dlist, brick.mims,
1866 brick.region, version);
1869 brick.pbr->complex_post_assembly_in_serial(*
this, ib, brick.vlist,
1870 brick.dlist, brick.mims,
1872 brick.cveclist[rhs_ind],
1873 brick.cveclist_sym[rhs_ind],
1874 brick.region, version);
1876 if (brick.is_update_brick)
1878 for (
auto &&mat : brick.cmatlist)
1881 for (
auto &&vecs : brick.cveclist)
1882 for (
auto &&vec : vecs)
1885 for (
auto &&vecs : brick.cveclist_sym)
1886 for (
auto &&vec : vecs)
1892 brick.pbr->real_pre_assembly_in_serial(*
this, ib, brick.vlist,
1893 brick.dlist, brick.mims,
1895 brick.rveclist[rhs_ind],
1896 brick.rveclist_sym[rhs_ind],
1897 brick.region, version);
1900 accumulated_distro<real_matlist> rmatlist(brick.rmatlist);
1901 accumulated_distro<real_veclist> rveclist(brick.rveclist[rhs_ind]);
1902 accumulated_distro<real_veclist> rveclist_sym(brick.rveclist_sym[rhs_ind]);
1906 brick.pbr->asm_real_tangent_terms(*
this, ib, brick.vlist,
1907 brick.dlist, brick.mims,
1915 brick.pbr->real_post_assembly_in_serial(*
this, ib, brick.vlist,
1916 brick.dlist, brick.mims,
1918 brick.rveclist[rhs_ind],
1919 brick.rveclist_sym[rhs_ind],
1920 brick.region, version);
1922 if (brick.is_update_brick)
1924 for (
auto &&mat : brick.rmatlist)
1927 for (
auto &&vecs : brick.rveclist)
1928 for (
auto &&vec : vecs)
1931 for (
auto &&vecs : brick.rveclist_sym)
1932 for (
auto &&vec : vecs)
1939 void model::set_dispatch_coeff() {
1940 for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) {
1941 brick_description &brick = bricks[ib];
1942 if (brick.pdispatch)
1943 brick.pdispatch->set_dispatch_coeff(*
this, ib);
1949 context_check();
if (act_size_to_be_done) actualize_sizes();
1950 for (
auto && v : variables) v.second.clear_temporaries();
1952 set_dispatch_coeff();
1954 for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) {
1955 brick_description &brick = bricks[ib];
1956 if (brick.pdispatch) {
1958 brick.pdispatch->next_complex_iter(*
this, ib, brick.vlist,
1960 brick.cmatlist, brick.cveclist,
1961 brick.cveclist_sym,
true);
1963 brick.pdispatch->next_real_iter(*
this, ib, brick.vlist, brick.dlist,
1964 brick.rmatlist, brick.rveclist,
1965 brick.rveclist_sym,
true);
1971 context_check();
if (act_size_to_be_done) actualize_sizes();
1972 set_dispatch_coeff();
1974 for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) {
1975 brick_description &brick = bricks[ib];
1976 if (brick.pdispatch) {
1978 brick.pdispatch->next_complex_iter(*
this, ib, brick.vlist,
1980 brick.cmatlist, brick.cveclist,
1981 brick.cveclist_sym,
false);
1983 brick.pdispatch->next_real_iter(*
this, ib, brick.vlist, brick.dlist,
1984 brick.rmatlist, brick.rveclist,
1985 brick.rveclist_sym,
false);
1989 for (
auto &&v : variables)
1990 for (
size_type i = 1; i < v.second.n_iter; ++i) {
1992 gmm::copy(v.second.complex_value[i-1], v.second.complex_value[i]);
1994 gmm::copy(v.second.real_value[i-1], v.second.real_value[i]);
1995 v.second.v_num_data[i] = act_counter();
1999 bool model::is_var_newer_than_brick(
const std::string &varname,
2001 const brick_description &brick = bricks[ib];
2002 var_description &vd = variables[varname];
2003 if (niter ==
size_type(-1)) niter = vd.default_iter;
2004 return (vd.v_num > brick.v_num || vd.v_num_data[niter] > brick.v_num);
2007 bool model::is_var_mf_newer_than_brick(
const std::string &varname,
2009 const brick_description &brick = bricks[ib];
2010 var_description &vd = variables[varname];
2011 return (vd.v_num > brick.v_num);
2014 bool model::is_mim_newer_than_brick(
const mesh_im &im,
2016 const brick_description &brick = bricks[ib];
2017 return (im.version_number() > brick.v_num);
2020 void model::define_variable_group(
const std::string &group_name,
2021 const std::vector<std::string> &nl) {
2023 "variables cannot be the same as a variable name");
2025 std::set<const mesh *> ms;
2026 bool is_data_ =
false;
2027 for (
size_type i = 0; i < nl.size(); ++i) {
2032 "It is not possible to mix variables and data in a group");
2035 "All variables in a group have to exist in the model");
2037 GMM_ASSERT1(mf,
"Variables in a group should be fem variables");
2038 GMM_ASSERT1(ms.find(&(mf->linked_mesh())) == ms.end(),
2039 "Two variables in a group cannot share the same mesh");
2040 ms.insert(&(mf->linked_mesh()));
2042 variable_groups[group_name] = nl;
2045 void model::add_assembly_assignments(
const std::string &varname,
2048 GMM_ASSERT1(order < 3 || order ==
size_type(-1),
"Bad order value");
2049 const im_data *imd = pim_data_of_variable(varname);
2050 GMM_ASSERT1(imd != 0,
"Only applicable to im_data");
2051 assignement_desc as;
2052 as.varname = varname; as.expr = expr; as.region = rg; as.order = order;
2054 assignments.push_back(as);
2057 void model::add_temporaries(
const varnamelist &vl,
2058 gmm::uint64_type id_num)
const {
2059 for (
size_type i = 0; i < vl.size(); ++i) {
2060 var_description &vd = variables[vl[i]];
2061 if (vd.n_iter > 1) {
2062 vd.add_temporary(id_num);
2067 bool model::temporary_uptodate(
const std::string &varname,
2068 gmm::uint64_type id_num,
2070 var_description &vd = variables[varname];
2072 for (; ind < vd.n_iter + vd.n_temp_iter ; ++ind) {
2073 if (vd.v_num_var_iter[ind] == id_num)
break;
2075 if (ind < vd.n_iter + vd.n_temp_iter) {
2076 if (vd.v_num_iter[ind] <= vd.v_num_data[vd.default_iter]) {
2077 vd.v_num_iter[ind] = act_counter();
2086 void model::set_default_iter_of_variable(
const std::string &varname,
2089 var_description &vd = variables[varname];
2090 GMM_ASSERT1(ind < vd.n_iter + vd.n_temp_iter,
2091 "Inexistent iteration " << ind);
2092 vd.default_iter = ind;
2096 void model::reset_default_iter_of_variables(
const varnamelist &vl)
const {
2097 for (
size_type i = 0; i < vl.size(); ++i)
2098 variables[vl[i]].default_iter = 0;
2101 const model_real_sparse_matrix &
2103 GMM_ASSERT1(bricks[ib].tlist[iterm].is_matrix_term,
2104 "Not a matrix term !");
2105 GMM_ASSERT1(bricks[ib].pbr->is_linear(),
"Nonlinear term !");
2106 return bricks[ib].rmatlist[iterm];
2109 const model_complex_sparse_matrix &
2111 GMM_ASSERT1(bricks[ib].tlist[iterm].is_matrix_term,
2112 "Not a matrix term !");
2113 GMM_ASSERT1(bricks[ib].pbr->is_linear(),
"Nonlinear term !");
2114 return bricks[ib].cmatlist[iterm];
2118 void model::update_brick(
size_type ib, build_version version)
const {
2119 const brick_description &brick = bricks[ib];
2120 bool cplx =
is_complex() && brick.pbr->is_complex();
2121 bool tobecomputed = brick.terms_to_be_computed
2122 || brick.pbr->is_to_be_computed_each_time()
2123 || !(brick.pbr->is_linear());
2126 if (!tobecomputed ) {
2127 for (
size_type i = 0; i < brick.vlist.size() && !tobecomputed; ++i) {
2128 var_description &vd = variables[brick.vlist[i]];
2129 if (vd.v_num > brick.v_num) tobecomputed =
true;
2134 for (
size_type i = 0; i < brick.dlist.size() && !tobecomputed; ++i) {
2135 var_description &vd = variables[brick.dlist[i]];
2136 if (vd.v_num > brick.v_num || vd.v_num_data[vd.default_iter] > brick.v_num) {
2137 tobecomputed =
true;
2138 version = build_version(version | BUILD_ON_DATA_CHANGE);
2143 if (!tobecomputed ) {
2144 for (
size_type i = 0; i < brick.mims.size() && !tobecomputed; ++i) {
2145 if (brick.mims[i]->version_number() > brick.v_num) tobecomputed =
true;
2150 brick.external_load = scalar_type(0);
2152 if (!(brick.pdispatch))
2153 { brick_call(ib, version, 0); }
2156 brick.pdispatch->asm_complex_tangent_terms
2157 (*
this, ib, brick.cmatlist, brick.cveclist, brick.cveclist_sym,
2160 brick.pdispatch->asm_real_tangent_terms
2161 (*
this, ib, brick.rmatlist, brick.rveclist, brick.rveclist_sym,
2164 brick.v_num = act_counter();
2167 if (brick.pbr->is_linear()) brick.terms_to_be_computed =
false;
2174 const brick_description &brick = bricks[ib];
2175 if (brick.pbr->is_linear()) {
2177 bool cplx =
is_complex() && brick.pbr->is_complex();
2179 for (
size_type j = 0; j < brick.tlist.size(); ++j) {
2180 const term_description &term = brick.tlist[j];
2181 bool isg = term.is_global;
2184 size_type n_iter_1 = n_iter, n_iter_2 = n_iter;
2186 n_iter_1 = variables[term.var1].default_iter;
2187 if (term.is_matrix_term)
2188 n_iter_2 = variables[term.var2].default_iter;
2193 if (term.is_matrix_term) {
2196 model_complex_plain_vector V(nbgdof);
2197 for (VAR_SET::iterator it = variables.begin();
2198 it != variables.end(); ++it)
2199 if (it->second.is_variable) {
2201 ? it->second.default_iter : n_iter;
2202 gmm::copy(it->second.complex_value[n_iter_i],
2203 gmm::sub_vector(V, it->second.I));
2207 gmm::scaled(V, complex_type(-1)),
2208 brick.cveclist[ind_data][j]);
2212 gmm::scaled(variables[term.var2].complex_value[n_iter_2],
2214 brick.cveclist[ind_data][j]);
2218 model_real_plain_vector V(nbgdof);
2219 for (VAR_SET::iterator it = variables.begin();
2220 it != variables.end(); ++it)
2221 if (it->second.is_variable) {
2223 ? it->second.default_iter : n_iter;
2224 gmm::copy(it->second.real_value[n_iter_i],
2225 gmm::sub_vector(V, it->second.I));
2228 (brick.rmatlist[j], gmm::scaled(V, scalar_type(-1)),
2229 brick.rveclist[ind_data][j]);
2233 gmm::scaled(variables[term.var2].real_value[n_iter_2],
2234 scalar_type(-1)), brick.rveclist[ind_data][j]);
2237 if (term.is_symmetric && term.var1.compare(term.var2)) {
2240 (gmm::conjugated(brick.cmatlist[j]),
2241 gmm::scaled(variables[term.var1].complex_value[n_iter_1],
2243 brick.cveclist_sym[ind_data][j]);
2246 (gmm::transposed(brick.rmatlist[j]),
2247 gmm::scaled(variables[term.var1].real_value[n_iter_1],
2249 brick.rveclist_sym[ind_data][j]);
2256 void model::update_affine_dependent_variables() {
2257 for (VAR_SET::iterator it = variables.begin(); it != variables.end(); ++it)
2258 if (it->second.is_affine_dependent) {
2259 VAR_SET::iterator it2 = variables.find(it->second.org_name);
2260 if (it->second.size() != it2->second.size())
2261 it->second.set_size();
2262 if (it->second.is_complex) {
2263 gmm::add(gmm::scaled(it2->second.complex_value[0],
2264 complex_type(it->second.alpha)),
2265 it->second.affine_complex_value,
2266 it->second.complex_value[0]);
2268 gmm::add(gmm::scaled(it2->second.real_value[0], it->second.alpha),
2269 it->second.affine_real_value, it->second.real_value[0]);
2271 it->second.v_num = std::max(it->second.v_num, it2->second.v_num);
2272 for (
size_type i = 0; i < it->second.n_iter; ++i)
2274 it->second.v_num_data[i] = std::max(it->second.v_num_data[i],
2275 it2->second.v_num_data[i]);
2286 GMM_ASSERT1(mf,
"Works only with fem variables.");
2290 for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) {
2291 brick_description &brick = bricks[ib];
2293 bool detected =
false;
2294 for (
size_type i = 0; i < brick.vlist.size(); ++i)
2295 if (brick.vlist[i].compare(varname) == 0)
2296 { detected =
true;
break; }
2298 if (detected && brick.mims.size()) {
2300 for (
auto &pmim : brick.mims)
2303 pmim->linked_mesh()));
2304 GMM_ASSERT1(ifo >= 0,
2305 "The given region is only partially covered by "
2306 "region of brick \"" << brick.pbr->brick_name()
2307 <<
"\". Please subdivise the region");
2309 std::string expr = brick.pbr->declare_volume_assembly_string
2310 (*
this, ib, brick.vlist, brick.dlist);
2312 ga_workspace workspace(*
this);
2313 size_type order = workspace.add_expression
2314 (expr, dummy_mim, region);
2315 GMM_ASSERT1(order <= 1,
"Wrong order for a Neumann term");
2316 expr = workspace.extract_Neumann_term(varname);
2319 result +=
" + " + workspace.extract_Neumann_term(varname);
2321 result = workspace.extract_Neumann_term(varname);
2333 GMM_ASSERT1(version != BUILD_ON_DATA_CHANGE,
2334 "Invalid assembly version BUILD_ON_DATA_CHANGE");
2335 GMM_ASSERT1(version != BUILD_WITH_LIN,
2336 "Invalid assembly version BUILD_WITH_LIN");
2337 GMM_ASSERT1(version != BUILD_WITH_INTERNAL,
2338 "Invalid assembly version BUILD_WITH_INTERNAL");
2340 #if GETFEM_PARA_LEVEL > 0
2341 double t_ref = MPI_Wtime();
2343 MPI_Comm_rank(MPI_COMM_WORLD, &rk);
2344 MPI_Comm_size(MPI_COMM_WORLD, &nbp);
2347 context_check();
if (act_size_to_be_done) actualize_sizes();
2349 if (version & BUILD_MATRIX) gmm::clear(cTM);
2350 if (version & BUILD_RHS) gmm::clear(crhs);
2353 if (version & BUILD_MATRIX) gmm::clear(rTM);
2354 if (version & BUILD_RHS) gmm::clear(rrhs);
2356 clear_dof_constraints();
2357 generic_expressions.clear();
2358 update_affine_dependent_variables();
2360 if (version & BUILD_RHS) approx_external_load_ = scalar_type(0);
2362 for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) {
2364 brick_description &brick = bricks[ib];
2367 bool auto_disabled_brick =
true;
2368 for (
size_type j = 0; j < brick.vlist.size(); ++j) {
2370 auto_disabled_brick =
false;
2372 if (auto_disabled_brick)
continue;
2374 update_brick(ib, version);
2376 bool cplx =
is_complex() && brick.pbr->is_complex();
2378 scalar_type coeff0 = scalar_type(1);
2379 if (brick.pdispatch) coeff0 = brick.matrix_coeff;
2383 for (
size_type j = 0; j < brick.tlist.size(); ++j) {
2384 term_description &term = brick.tlist[j];
2385 bool isg = term.is_global, isprevious =
false;
2387 scalar_type alpha = coeff0, alpha1 = coeff0, alpha2 = coeff0;
2388 gmm::sub_interval I1(0,nbgdof), I2(0,nbgdof);
2389 var_description *var1=
nullptr, *var2=
nullptr;
2391 VAR_SET::iterator it1 = variables.find(term.var1);
2392 GMM_ASSERT1(it1 != variables.end(),
"Internal error");
2393 var1 = &(it1->second);
2394 GMM_ASSERT1(var1->is_variable,
"Assembly of data not allowed");
2396 if (term.is_matrix_term) {
2397 VAR_SET::iterator it2 = variables.find(term.var2);
2398 GMM_ASSERT1(it2 != variables.end(),
"Internal error");
2399 var2 = &(it2->second);
2401 if (!(var2->is_variable)) {
2402 std::string vorgname = sup_previous_and_dot_to_varname(term.var2);
2403 VAR_SET::iterator it3 = variables.find(vorgname);
2404 GMM_ASSERT1(it3->second.is_variable,
2405 "Assembly of data not allowed");
2409 alpha *= var1->alpha * var2->alpha;
2410 alpha1 *= var1->alpha;
2411 alpha2 *= var2->alpha;
2416 if (term.is_matrix_term && (version & BUILD_MATRIX) && !isprevious
2417 && (isg || (var1->is_enabled() && var2->is_enabled()))) {
2418 gmm::add(gmm::scaled(brick.cmatlist[j], alpha),
2419 gmm::sub_matrix(cTM, I1, I2));
2420 if (term.is_symmetric && I1.first() != I2.first())
2421 gmm::add(gmm::scaled(gmm::transposed(brick.cmatlist[j]), alpha),
2422 gmm::sub_matrix(cTM, I2, I1));
2424 if (version & BUILD_RHS) {
2426 if (isg || var1->is_enabled()) {
2427 if (brick.pdispatch)
2428 for (
size_type k = 0; k < brick.nbrhs; ++k)
2429 gmm::add(gmm::scaled(brick.cveclist[k][j],
2431 gmm::sub_vector(crhs, I1));
2433 gmm::add(gmm::scaled(brick.cveclist[0][j],
2434 complex_type(alpha1)),
2435 gmm::sub_vector(crhs, I1));
2437 if (term.is_matrix_term && brick.pbr->is_linear() &&
is_linear()) {
2438 if (var2->is_affine_dependent && var1->is_enabled())
2439 gmm::mult_add(brick.cmatlist[j],
2440 gmm::scaled(var2->affine_complex_value,
2441 complex_type(-alpha1)),
2442 gmm::sub_vector(crhs, I1));
2443 if (term.is_symmetric && I1.first() != I2.first()
2444 && var1->is_affine_dependent && var2->is_enabled())
2445 gmm::mult_add(gmm::conjugated(brick.cmatlist[j]),
2446 gmm::scaled(var1->affine_complex_value,
2447 complex_type(-alpha2)),
2448 gmm::sub_vector(crhs, I2));
2450 if (term.is_matrix_term && brick.pbr->is_linear()
2451 && (!
is_linear() || (version & BUILD_WITH_LIN))
2452 && var1->is_enabled())
2453 gmm::mult_add(brick.cmatlist[j],
2454 gmm::scaled(var2->complex_value[0],
2455 complex_type(-alpha1)),
2456 gmm::sub_vector(crhs, I1));
2457 if (term.is_symmetric && I1.first() != I2.first()
2458 && var2->is_enabled()) {
2459 if (brick.pdispatch)
2460 for (
size_type k = 0; k < brick.nbrhs; ++k)
2461 gmm::add(gmm::scaled(brick.cveclist_sym[k][j],
2463 gmm::sub_vector(crhs, I2));
2465 gmm::add(gmm::scaled(brick.cveclist_sym[0][j],
2466 complex_type(alpha2)),
2467 gmm::sub_vector(crhs, I2));
2468 if (brick.pbr->is_linear()
2469 && (!
is_linear() || (version & BUILD_WITH_LIN)))
2470 gmm::mult_add(gmm::conjugated(brick.cmatlist[j]),
2471 gmm::scaled(var1->complex_value[0],
2472 complex_type(-alpha2)),
2473 gmm::sub_vector(crhs, I2));
2477 if (term.is_matrix_term && (version & BUILD_MATRIX) && !isprevious
2478 && (isg || (var1->is_enabled() && var2->is_enabled()))) {
2479 gmm::add(gmm::scaled(brick.rmatlist[j], alpha),
2480 gmm::sub_matrix(cTM, I1, I2));
2481 if (term.is_symmetric && I1.first() != I2.first())
2482 gmm::add(gmm::scaled(gmm::transposed(brick.rmatlist[j]), alpha),
2483 gmm::sub_matrix(cTM, I2, I1));
2485 if (version & BUILD_RHS) {
2487 if (isg || var1->is_enabled()) {
2488 if (brick.pdispatch)
2489 for (
size_type k = 0; k < brick.nbrhs; ++k)
2490 gmm::add(gmm::scaled(brick.rveclist[k][j],
2492 gmm::sub_vector(crhs, I1));
2494 gmm::add(gmm::scaled(brick.rveclist[0][j], alpha1),
2495 gmm::sub_vector(crhs, I1));
2497 if (term.is_matrix_term && brick.pbr->is_linear() &&
is_linear()) {
2498 if (var2->is_affine_dependent && var1->is_enabled())
2499 gmm::mult_add(brick.rmatlist[j],
2500 gmm::scaled(var2->affine_complex_value,
2501 complex_type(-alpha1)),
2502 gmm::sub_vector(crhs, I1));
2503 if (term.is_symmetric && I1.first() != I2.first()
2504 && var1->is_affine_dependent && var2->is_enabled())
2505 gmm::mult_add(gmm::transposed(brick.rmatlist[j]),
2506 gmm::scaled(var1->affine_complex_value,
2507 complex_type(-alpha2)),
2508 gmm::sub_vector(crhs, I2));
2510 if (term.is_matrix_term && brick.pbr->is_linear()
2511 && (!
is_linear() || (version & BUILD_WITH_LIN))
2512 && var1->is_enabled())
2513 gmm::mult_add(brick.rmatlist[j],
2514 gmm::scaled(var2->complex_value[0],
2515 complex_type(-alpha1)),
2516 gmm::sub_vector(crhs, I1));
2517 if (term.is_symmetric && I1.first() != I2.first()
2518 && var2->is_enabled()) {
2519 if (brick.pdispatch)
2520 for (
size_type k = 0; k < brick.nbrhs; ++k)
2521 gmm::add(gmm::scaled(brick.rveclist_sym[k][j],
2523 gmm::sub_vector(crhs, I2));
2525 gmm::add(gmm::scaled(brick.rveclist_sym[0][j], alpha2),
2526 gmm::sub_vector(crhs, I2));
2528 if (brick.pbr->is_linear()
2529 && (!
is_linear() || (version & BUILD_WITH_LIN)))
2530 gmm::mult_add(gmm::transposed(brick.rmatlist[j]),
2531 gmm::scaled(var1->complex_value[0],
2532 complex_type(-alpha2)),
2533 gmm::sub_vector(crhs, I2));
2537 if (term.is_matrix_term && (version & BUILD_MATRIX) && !isprevious
2538 && (isg || (var1->is_enabled() && var2->is_enabled()))) {
2539 gmm::add(gmm::scaled(brick.rmatlist[j], alpha),
2540 gmm::sub_matrix(rTM, I1, I2));
2541 if (term.is_symmetric && I1.first() != I2.first())
2542 gmm::add(gmm::scaled(gmm::transposed(brick.rmatlist[j]), alpha),
2543 gmm::sub_matrix(rTM, I2, I1));
2545 if (version & BUILD_RHS) {
2547 auto vec_out1 = gmm::sub_vector(rrhs, I1);
2548 if (isg || var1->is_enabled()) {
2549 if (brick.pdispatch)
2550 for (
size_type k = 0; k < brick.nbrhs; ++k)
2551 gmm::add(gmm::scaled(brick.rveclist[k][j],
2555 gmm::add(gmm::scaled(brick.rveclist[0][j], alpha1),
2558 if (var1->is_enabled()
2559 && term.is_matrix_term && brick.pbr->is_linear()) {
2560 bool affine_contrib(
is_linear() && var2->is_affine_dependent);
2561 bool linear_contrib(!
is_linear() || (version & BUILD_WITH_LIN));
2562 const auto &matj = brick.rmatlist[j];
2563 const auto vec_affine2 = gmm::scaled(var2->affine_real_value,
2565 const auto vec_linear2 = gmm::scaled(var2->real_value[0],
2568 model_real_plain_vector vec_tmp1(I1.size(), 0.);
2570 gmm::mult(matj, vec_affine2, vec_tmp1);
2572 gmm::mult_add(matj, vec_linear2, vec_tmp1);
2573 MPI_SUM_VECTOR(vec_tmp1);
2574 gmm::add(vec_tmp1, vec_out1);
2577 gmm::mult_add(matj, vec_affine2, vec_out1);
2579 gmm::mult_add(matj, vec_linear2, vec_out1);
2584 if (term.is_symmetric && I1.first() != I2.first() &&
2585 var2->is_enabled()) {
2586 auto vec_out2 = gmm::sub_vector(rrhs, I2);
2587 if (brick.pdispatch)
2588 for (
size_type k = 0; k < brick.nbrhs; ++k)
2589 gmm::add(gmm::scaled(brick.rveclist_sym[k][j],
2593 gmm::add(gmm::scaled(brick.rveclist_sym[0][j], alpha2),
2595 if (term.is_matrix_term && brick.pbr->is_linear()) {
2596 bool affine_contrib(
is_linear() && var1->is_affine_dependent);
2597 bool linear_contrib(!
is_linear() || (version & BUILD_WITH_LIN));
2598 const auto matj_trans = gmm::transposed(brick.rmatlist[j]);
2599 const auto vec_affine1 = gmm::scaled(var1->affine_real_value,
2601 const auto vec_linear1 = gmm::scaled(var1->real_value[0],
2604 model_real_plain_vector vec_tmp2(I2.size(),0.);
2606 gmm::mult(matj_trans, vec_affine1, vec_tmp2);
2608 gmm::mult_add(matj_trans, vec_linear1, vec_tmp2);
2609 MPI_SUM_VECTOR(vec_tmp2);
2610 gmm::add(vec_tmp2, vec_out2);
2613 gmm::mult_add(matj_trans, vec_affine1, vec_out2);
2615 gmm::mult_add(matj_trans, vec_linear1, vec_out2);
2623 if (brick.pbr->is_linear())
2624 brick.terms_to_be_computed =
false;
2636 if (version & BUILD_RHS) approx_external_load_ += brick.external_load;
2639 if (version & BUILD_RHS && version & BUILD_WITH_INTERNAL) {
2642 gmm::sub_interval IP(0,gmm::vect_size(rrhs));
2643 gmm::fill(full_rrhs, 0.);
2644 gmm::copy(rrhs, gmm::sub_vector(full_rrhs, IP));
2648 if (generic_expressions.size()) {
2651 if (version & BUILD_RHS)
2652 GMM_TRACE2(
"Global generic assembly RHS");
2653 if (version & BUILD_MATRIX)
2654 GMM_TRACE2(
"Global generic assembly tangent term");
2657 auto add_assignments_and_expressions_to_workspace =
2658 [&](ga_workspace &workspace)
2660 for (
const auto &ad : assignments)
2661 workspace.add_assignment_expression
2662 (ad.varname, ad.expr, ad.region, ad.order, ad.before);
2663 for (
const auto &ge : generic_expressions)
2664 workspace.add_expression(ge.expr, ge.mim, ge.region,
2665 2, ge.secondary_domain);
2668 const bool with_internal = version & BUILD_WITH_INTERNAL
2670 model_real_sparse_matrix intern_mat;
2671 model_real_plain_vector res0,
2674 size_type full_size = gmm::vect_size(full_rrhs),
2675 primary_size = gmm::vect_size(rrhs);
2677 if ((version & BUILD_RHS) || (version & BUILD_MATRIX && with_internal))
2678 gmm::resize(res0, with_internal ? full_size : primary_size);
2679 if (version & BUILD_MATRIX && with_internal)
2680 gmm::resize(res1, full_size);
2682 if (version & BUILD_MATRIX) {
2683 if (with_internal) {
2684 gmm::resize(intern_mat, full_size, primary_size);
2685 gmm::resize(res1, full_size);
2691 if (version & BUILD_RHS) {
2694 ga_workspace workspace(*
this);
2695 add_assignments_and_expressions_to_workspace(workspace);
2696 workspace.set_assembled_vector(res0_distro);
2697 workspace.assembly(1, with_internal);
2698 if (with_internal) {
2699 gmm::copy(res0_distro.get(), res1_distro.get());
2700 gmm::add(gmm::scaled(full_rrhs, scalar_type(-1)),
2702 workspace.set_assembled_vector(res1_distro);
2703 workspace.set_internal_coupling_matrix(intern_mat_distro);
2705 workspace.set_assembled_matrix(tangent_matrix_distro);
2706 workspace.assembly(2, with_internal);
2711 ga_workspace workspace(*
this);
2712 add_assignments_and_expressions_to_workspace(workspace);
2713 if (with_internal) {
2714 gmm::copy(gmm::scaled(full_rrhs, scalar_type(-1)),
2716 workspace.set_assembled_vector(res1_distro);
2717 workspace.set_internal_coupling_matrix(intern_mat_distro);
2719 workspace.set_assembled_matrix(tangent_matrix_distro);
2720 workspace.assembly(2, with_internal);
2724 else if (version & BUILD_RHS) {
2727 ga_workspace workspace(*
this);
2728 add_assignments_and_expressions_to_workspace(workspace);
2729 workspace.set_assembled_vector(res0_distro);
2730 workspace.assembly(1, with_internal);
2734 if (version & BUILD_RHS) {
2735 gmm::scale(res0, scalar_type(-1));
2736 if (with_internal) {
2737 gmm::sub_interval IP(0,gmm::vect_size(rrhs));
2738 gmm::add(gmm::sub_vector(res0, IP), rrhs);
2739 gmm::add(res0, full_rrhs);
2741 gmm::add(res0, rrhs);
2744 if (version & BUILD_MATRIX && with_internal) {
2745 gmm::scale(res1, scalar_type(-1));
2746 gmm::sub_interval IP(0, primary_size),
2747 II(primary_size, full_size-primary_size);
2748 gmm::copy(gmm::sub_matrix(intern_mat, II, IP), internal_rTM);
2749 gmm::add(gmm::sub_vector(res1, IP), rrhs);
2750 gmm::copy(gmm::sub_vector(res1, II), internal_sol);
2755 if ((version & BUILD_RHS) || (version & BUILD_MATRIX)) {
2757 std::vector<size_type> dof_indices;
2758 std::vector<complex_type> dof_pr_values;
2759 std::vector<complex_type> dof_go_values;
2760 for (
const auto &keyval : complex_dof_constraints) {
2761 const gmm::sub_interval &I = interval_of_variable(keyval.first);
2763 for (
const auto &val : keyval.second) {
2764 dof_indices.push_back(val.first + I.first());
2765 dof_go_values.push_back(val.second);
2766 dof_pr_values.push_back(V[val.first]);
2770 if (dof_indices.size()) {
2771 gmm::sub_index SI(dof_indices);
2772 gmm::sub_interval II(0,
nb_dof());
2774 if (version & BUILD_RHS) {
2775 if (MPI_IS_MASTER())
2776 approx_external_load_ += gmm::vect_norm1(dof_go_values);
2778 if (is_symmetric_) {
2779 scalar_type valnorm = gmm::vect_norm2(dof_go_values);
2780 if (valnorm > scalar_type(0)) {
2781 GMM_ASSERT1(version & BUILD_MATRIX,
"Rhs only for a "
2782 "symmetric linear problem with dof "
2783 "constraint not allowed");
2784 model_complex_plain_vector vv(gmm::vect_size(crhs));
2785 gmm::mult(gmm::sub_matrix(cTM, II, SI), dof_go_values, vv);
2787 gmm::add(gmm::scaled(vv, scalar_type(-1)), crhs);
2790 gmm::copy(dof_go_values, gmm::sub_vector(crhs, SI));
2792 gmm::add(dof_go_values,
2793 gmm::scaled(dof_pr_values, complex_type(-1)),
2794 gmm::sub_vector(crhs, SI));
2797 if (version & BUILD_MATRIX) {
2798 gmm::clear(gmm::sub_matrix(cTM, SI, II));
2799 if (is_symmetric_) gmm::clear(gmm::sub_matrix(cTM, II, SI));
2801 if (MPI_IS_MASTER()) {
2802 for (
size_type i = 0; i < dof_indices.size(); ++i)
2803 cTM(dof_indices[i], dof_indices[i]) = complex_type(1);
2808 std::vector<size_type> dof_indices;
2809 std::vector<scalar_type> dof_pr_values;
2810 std::vector<scalar_type> dof_go_values;
2811 for (
const auto &keyval : real_dof_constraints) {
2812 const gmm::sub_interval &I = interval_of_variable(keyval.first);
2813 const model_real_plain_vector &V =
real_variable(keyval.first);
2814 for (
const auto &val : keyval.second) {
2815 dof_indices.push_back(val.first + I.first());
2816 dof_go_values.push_back(val.second);
2817 dof_pr_values.push_back(V[val.first]);
2821 #if GETFEM_PARA_LEVEL > 1
2822 GMM_ASSERT1(MPI_IS_MASTER() || (dof_indices.size() == 0),
2823 "Sorry, for the moment, the dof constraints have to be "
2824 "added on the master process only");
2825 size_type dof_indices_size = dof_indices.size();
2826 MPI_BCAST0_SCALAR(dof_indices_size);
2827 dof_indices.resize(dof_indices_size);
2828 MPI_BCAST0_VECTOR(dof_indices);
2829 dof_pr_values.resize(dof_indices_size);
2830 MPI_BCAST0_VECTOR(dof_pr_values);
2831 dof_go_values.resize(dof_indices_size);
2832 MPI_BCAST0_VECTOR(dof_go_values);
2835 if (dof_indices.size()) {
2836 gmm::sub_index SI(dof_indices);
2837 gmm::sub_interval II(0,
nb_dof());
2839 if (version & BUILD_RHS) {
2840 if (MPI_IS_MASTER())
2841 approx_external_load_ += gmm::vect_norm1(dof_go_values);
2843 if (is_symmetric_) {
2844 scalar_type valnorm = gmm::vect_norm2(dof_go_values);
2845 if (valnorm > scalar_type(0)) {
2846 GMM_ASSERT1(version & BUILD_MATRIX,
"Rhs only for a "
2847 "symmetric linear problem with dof "
2848 "constraint not allowed");
2849 model_real_plain_vector vv(gmm::vect_size(rrhs));
2850 gmm::mult(gmm::sub_matrix(rTM, II, SI), dof_go_values, vv);
2852 gmm::add(gmm::scaled(vv, scalar_type(-1)), rrhs);
2855 gmm::copy(dof_go_values, gmm::sub_vector(rrhs, SI));
2857 gmm::add(dof_go_values,
2858 gmm::scaled(dof_pr_values, scalar_type(-1)),
2859 gmm::sub_vector(rrhs, SI));
2862 if (version & BUILD_MATRIX) {
2863 gmm::clear(gmm::sub_matrix(rTM, SI, II));
2864 if (is_symmetric_) gmm::clear(gmm::sub_matrix(rTM, II, SI));
2866 if (MPI_IS_MASTER()) {
2867 for (
size_type i = 0; i < dof_indices.size(); ++i)
2868 rTM(dof_indices[i], dof_indices[i]) = scalar_type(1);
2875 if (version & BUILD_RHS) {
2878 MPI_BCAST0_SCALAR(approx_external_load_);
2881 #if GETFEM_PARA_LEVEL > 0
2882 if (MPI_IS_MASTER()) cout <<
"Assembly time " << MPI_Wtime()-t_ref << endl;
2891 return it->second.associated_mf();
2897 return it->second.passociated_mf();
2901 model::qdims_of_variable(
const std::string &name)
const {
2903 const mesh_fem *mf = it->second.passociated_mf();
2904 const im_data *imd = it->second.imd;
2907 bgeot::multi_index mi = mf->get_qdims();
2908 if (n > 1 || it->second.qdims.size() > 1) {
2910 if (mi.back() == 1) { mi.back() *= it->second.qdims[0]; ++i; }
2911 for (; i < it->second.qdims.size(); ++i)
2912 mi.push_back(it->second.qdims[i]);
2916 bgeot::multi_index mi = imd->tensor_size();
2917 if (n > 1 || it->second.qdims.size() > 1) {
2919 if (mi.back() == 1) { mi.back() *= it->second.qdims[0]; ++i; }
2920 for (; i < it->second.qdims.size(); ++i)
2921 mi.push_back(it->second.qdims[i]);
2925 return it->second.qdims;
2928 size_type model::qdim_of_variable(
const std::string &name)
const {
2930 const mesh_fem *mf = it->second.passociated_mf();
2931 const im_data *imd = it->second.imd;
2934 return mf->get_qdim() * n;
2936 return imd->tensor_size().total_size() * n;
2942 const gmm::sub_interval &
2943 model::interval_of_variable(
const std::string &name)
const {
2945 if (act_size_to_be_done) actualize_sizes();
2946 VAR_SET::const_iterator it = find_variable(name);
2947 return it->second.I;
2950 const model_real_plain_vector &
2956 const model_real_plain_vector &
2958 GMM_ASSERT1(!complex_version,
"This model is a complex one");
2959 GMM_ASSERT1(!
is_old(name),
"Please don't use Old_ prefix in combination "
2960 "with variable version");
2962 auto it = variables.find(name);
2963 GMM_ASSERT1(it != variables.end(),
"Undefined variable " << name);
2964 if (act_size_to_be_done && it->second.mf) {
2965 if (it->second.filter != VDESCRFILTER_NO)
2968 it->second.set_size();
2970 if (niter ==
size_type(-1)) niter = it->second.default_iter;
2971 GMM_ASSERT1(it->second.n_iter + it->second.n_temp_iter > niter,
2972 "Invalid iteration number " << niter <<
" for " << name);
2973 return it->second.real_value[niter];
2976 const model_complex_plain_vector &
2982 const model_complex_plain_vector &
2984 GMM_ASSERT1(complex_version,
"This model is a real one");
2985 GMM_ASSERT1(!
is_old(name),
"Please don't use Old_ prefix in combination with"
2986 " variable version");
2988 auto it = variables.find(name);
2989 GMM_ASSERT1(it!=variables.end(),
"Undefined variable " << name);
2990 if (act_size_to_be_done && it->second.mf) {
2991 if (it->second.filter != VDESCRFILTER_NO)
2994 it->second.set_size();
2996 if (niter ==
size_type(-1)) niter = it->second.default_iter;
2997 GMM_ASSERT1(it->second.n_iter + it->second.n_temp_iter > niter,
2998 "Invalid iteration number "
2999 << niter <<
" for " << name);
3000 return it->second.complex_value[niter];
3003 model_real_plain_vector &
3010 model_real_plain_vector &
3012 GMM_ASSERT1(!complex_version,
"This model is a complex one");
3013 GMM_ASSERT1(!
is_old(name),
"Please don't use Old_ prefix in combination with"
3014 " variable version");
3016 auto it = variables.find(name);
3017 GMM_ASSERT1(it!=variables.end(),
"Undefined variable " << name);
3018 if (act_size_to_be_done && it->second.mf) {
3019 if (it->second.filter != VDESCRFILTER_NO)
3022 it->second.set_size();
3024 if (niter ==
size_type(-1)) niter = it->second.default_iter;
3025 it->second.v_num_data[niter] = act_counter();
3026 GMM_ASSERT1(it->second.n_iter + it->second.n_temp_iter > niter,
3027 "Invalid iteration number "
3028 << niter <<
" for " << name);
3029 return it->second.real_value[niter];
3032 model_complex_plain_vector &
3038 model_complex_plain_vector &
3040 GMM_ASSERT1(complex_version,
"This model is a real one");
3041 GMM_ASSERT1(!
is_old(name),
"Please don't use Old_ prefix in combination with"
3042 " variable version");
3044 auto it = variables.find(name);
3045 GMM_ASSERT1(it!=variables.end(),
"Undefined variable " << name);
3046 if (act_size_to_be_done && it->second.mf) {
3047 if (it->second.filter != VDESCRFILTER_NO)
3050 it->second.set_size();
3052 if (niter ==
size_type(-1)) niter = it->second.default_iter;
3053 it->second.v_num_data[niter] = act_counter();
3054 GMM_ASSERT1(it->second.n_iter + it->second.n_temp_iter > niter,
3055 "Invalid iteration number "
3056 << niter <<
" for " << name);
3057 return it->second.complex_value[niter];
3060 model_real_plain_vector &
3061 model::set_real_constant_part(
const std::string &name)
const {
3062 GMM_ASSERT1(!complex_version,
"This model is a complex one");
3064 VAR_SET::iterator it = variables.find(name);
3065 GMM_ASSERT1(it != variables.end(),
"Undefined variable " << name);
3066 GMM_ASSERT1(it->second.is_affine_dependent,
3067 "Only for affine dependent variables");
3068 if (act_size_to_be_done && it->second.mf) {
3069 if (it->second.filter != VDESCRFILTER_NO)
3072 it->second.set_size();
3074 for (
auto &v_num : it->second.v_num_data) v_num = act_counter();
3075 return it->second.affine_real_value;
3078 model_complex_plain_vector &
3079 model::set_complex_constant_part(
const std::string &name)
const {
3080 GMM_ASSERT1(complex_version,
"This model is a real one");
3082 VAR_SET::iterator it = variables.find(name);
3083 GMM_ASSERT1(it!=variables.end(),
"Undefined variable " << name);
3084 if (act_size_to_be_done && it->second.mf) {
3085 if (it->second.filter != VDESCRFILTER_NO)
3088 it->second.set_size();
3090 for (
auto &v_num : it->second.v_num_data) v_num = act_counter();
3091 return it->second.affine_complex_value;
3096 const brick_description &brick = bricks[ind_brick];
3097 update_brick(ind_brick, model::BUILD_ALL);
3099 brick.pbr->check_stiffness_matrix_and_rhs(*
this, ind_brick, brick.tlist,
3100 brick.vlist, brick.dlist, brick.mims, brick.rmatlist,
3101 brick.rveclist[0], brick.rveclist_sym[0], brick.region);
3104 void model::clear() {
3106 active_bricks.clear();
3107 valid_bricks.clear();
3108 real_dof_constraints.clear();
3109 complex_dof_constraints.clear();
3111 rTM = model_real_sparse_matrix();
3112 cTM = model_complex_sparse_matrix();
3113 rrhs = model_real_plain_vector();
3114 crhs = model_complex_plain_vector();
3119 void virtual_brick::full_asm_real_tangent_terms_(
const model &md,
size_type ind_brick,
3120 const model::varnamelist &term_list,
3121 const model::varnamelist &var_list,
3122 const model::mimlist &mim_list,
3123 model::real_matlist &rmatlist,
3124 model::real_veclist &rveclist,
3125 model::real_veclist &rveclist_sym,
3126 size_type region, build_version version)
const
3129 rveclist, rveclist_sym, region, version);
3131 rveclist, rveclist_sym, region, version);
3133 rveclist, rveclist_sym, region, version);
3138 const model::termlist& tlist,
3139 const model::varnamelist &vl,
3140 const model::varnamelist &dl,
3141 const model::mimlist &mims,
3142 model::real_matlist &matl,
3143 model::real_veclist &rvc1,
3144 model::real_veclist &rvc2,
3146 const scalar_type TINY)
const
3148 std::cout<<
"******Verifying stiffnesses of *******"<<std::endl;
3149 std::cout<<
"*** "<<brick_name()<<std::endl;
3152 std::map<std::string,size_type> rhs_index;
3153 for(
size_type iterm=0;iterm<matl.size();iterm++)
3154 if (tlist[iterm].var1==tlist[iterm].var2) rhs_index[tlist[iterm].var1]=iterm;
3156 if (rhs_index.size()==0){
3157 GMM_WARNING0(
"*** cannot verify stiffness for this brick***");
3160 full_asm_real_tangent_terms_(md, s, vl, dl, mims, matl, rvc1, rvc2,
3161 rg, model::BUILD_MATRIX);
3162 for(
size_type iterm=0;iterm<matl.size();iterm++)
3165 std::cout<<std::endl;
3166 std::cout<<
" Stiffness["<<tlist[iterm].var1
3167 <<
","<<tlist[iterm].var2<<
"]:"<<std::endl;
3170 std::cout<<
" "<<tlist[iterm].var1<<
" has zero size. Skipping this term"<<std::endl;
3175 std::cout<<
" "<<tlist[iterm].var2<<
" has zero size. Skipping this term"<<std::endl;
3179 model_real_sparse_matrix SM(matl[iterm]);
3180 gmm::fill(rvc1[rhs_index[tlist[iterm].var1]], 0.0);
3181 full_asm_real_tangent_terms_(md, s, vl, dl, mims, matl, rvc1, rvc2,
3182 rg, model::BUILD_RHS);
3183 if (gmm::mat_euclidean_norm(matl[iterm])<1e-12){
3184 std::cout<<
" The assembled matrix is nearly zero, skipping."<<std::endl;
3187 model_real_plain_vector RHS0(rvc1[rhs_index[tlist[iterm].var1]]);
3190 model_real_sparse_matrix fdSM(matl[iterm].nrows(), matl[iterm].ncols());
3192 model_real_plain_vector& RHS1 = rvc1[rhs_index[tlist[iterm].var1]];
3194 scalar_type relative_tiny = gmm::vect_norminf(RHS1)*TINY;
3195 if (relative_tiny < TINY) relative_tiny = TINY;
3197 for (
size_type j = 0; j < matl[iterm].ncols(); j++){
3198 U[j] += relative_tiny;
3199 gmm::fill(RHS1, 0.0);
3200 full_asm_real_tangent_terms_(md, s, vl, dl, mims, matl, rvc1, rvc2,
3201 rg, model::BUILD_RHS);
3202 for (
size_type i = 0; i<matl[iterm].nrows(); i++)
3203 fdSM(i, j) = (RHS0[i] - RHS1[i]) / relative_tiny;
3204 U[j] -= relative_tiny;
3207 model_real_sparse_matrix diffSM(matl[iterm].nrows(),matl[iterm].ncols());
3208 gmm::add(SM,gmm::scaled(fdSM,-1.0),diffSM);
3209 scalar_type norm_error_euc
3210 = gmm::mat_euclidean_norm(diffSM)/gmm::mat_euclidean_norm(SM)*100;
3211 scalar_type norm_error_1
3212 = gmm::mat_norm1(diffSM)/gmm::mat_norm1(SM)*100;
3213 scalar_type norm_error_max
3214 = gmm::mat_maxnorm(diffSM)/gmm::mat_maxnorm(SM)*100;
3217 scalar_type nsym_norm_error_euc=0.0;
3218 scalar_type nsym_norm_error_1=0.0;
3219 scalar_type nsym_norm_error_max=0.0;
3220 if (tlist[iterm].var1==tlist[iterm].var2){
3221 model_real_sparse_matrix diffSMtransposed(matl[iterm].nrows(),matl[iterm].ncols());
3222 gmm::add(gmm::transposed(fdSM),gmm::scaled(fdSM,-1.0),diffSMtransposed);
3224 = gmm::mat_euclidean_norm(diffSMtransposed)/gmm::mat_euclidean_norm(fdSM)*100;
3226 = gmm::mat_norm1(diffSMtransposed)/gmm::mat_norm1(fdSM)*100;
3228 = gmm::mat_maxnorm(diffSMtransposed)/gmm::mat_maxnorm(fdSM)*100;
3232 if(rvc1[0].size()<8){
3233 std::cout <<
"RHS Stiffness Matrix: \n";
3234 std::cout <<
"------------------------\n";
3235 for(
size_type i=0; i < rvc1[iterm].size(); ++i){
3237 for(
size_type j=0; j < rvc1[iterm].size(); ++j){
3238 std::cout << fdSM(i,j) <<
" ";
3242 std::cout <<
"Analytical Stiffness Matrix: \n";
3243 std::cout <<
"------------------------\n";
3244 for(
size_type i=0; i < rvc1[iterm].size(); ++i){
3246 for(
size_type j=0; j < rvc1[iterm].size(); ++j){
3247 std::cout << matl[iterm](i,j) <<
" ";
3251 std::cout <<
"Vector U: \n";
3252 std::cout <<
"------------------------\n";
3253 for(
size_type i=0; i < rvc1[iterm].size(); ++i){
3260 <<
"\n\nfinite diff test error_norm_eucl: " << norm_error_euc <<
"%\n"
3261 <<
"finite diff test error_norm1: " << norm_error_1 <<
"%\n"
3262 <<
"finite diff test error_max_norm: " << norm_error_max <<
"%\n\n\n";
3264 if (tlist[iterm].var1==tlist[iterm].var2){
3266 <<
"Nonsymmetrical test error_norm_eucl: "<< nsym_norm_error_euc<<
"%\n"
3267 <<
"Nonsymmetrical test error_norm1: " << nsym_norm_error_1 <<
"%\n"
3268 <<
"Nonsymmetrical test error_max_norm: " << nsym_norm_error_max <<
"%"
3288 struct gen_source_term_assembly_brick :
public virtual_brick {
3290 std::string expr, directvarname, directdataname;
3291 model::varnamelist vl_test1;
3292 std::string secondary_domain;
3295 const model::varnamelist &,
3296 const model::varnamelist &,
3297 const model::mimlist &mims,
3298 model::real_matlist &,
3299 model::real_veclist &vecl,
3300 model::real_veclist &,
3302 build_version)
const override {
3303 GMM_ASSERT1(vecl.size() == vl_test1.size()
3304 + ((directdataname.size() == 0) ? 0 : 1),
"Wrong number "
3305 "of terms for Generic source term assembly brick ");
3306 GMM_ASSERT1(mims.size() == 1,
"Generic source term assembly brick "
3307 "needs one and only one mesh_im");
3308 GMM_TRACE2(
"Generic source term assembly");
3310 gmm::clear(vecl[0]);
3313 const mesh_im &mim = *mims[0];
3318 ga_workspace workspace(md, ga_workspace::inherit::ALL);
3319 GMM_TRACE2(name <<
": generic source term assembly");
3320 workspace.add_expression(expr, mim, rg, 1, secondary_domain);
3321 workspace.assembly(1);
3322 const auto &V=workspace.assembled_vector();
3323 for (
size_type i = 0; i < vl_test1.size(); ++i) {
3324 const auto &I=workspace.interval_of_variable(vl_test1[i]);
3325 gmm::copy(gmm::sub_vector(V, I), vecl[i]);
3329 if (directvarname.size()) {
3334 void real_post_assembly_in_serial(
const model &md,
size_type ib,
3335 const model::varnamelist &,
3336 const model::varnamelist &,
3337 const model::mimlist &,
3338 model::real_matlist &,
3339 model::real_veclist &vecl,
3340 model::real_veclist &,
3342 build_version)
const override {
3343 scalar_type el = scalar_type(0);
3344 for (
size_type i=0; i < vecl.size(); ++i) el += gmm::vect_norm1(vecl[i]);
3345 md.add_external_load(ib, el);
3348 virtual std::string declare_volume_assembly_string
3349 (
const model &,
size_type,
const model::varnamelist &,
3350 const model::varnamelist &)
const {
3351 return std::string();
3354 gen_source_term_assembly_brick(
const std::string &expr_,
3355 std::string brickname,
3356 const model::varnamelist &vl_test1_,
3357 const std::string &directvarname_,
3358 const std::string &directdataname_,
3359 const std::string &secdom)
3360 : vl_test1(vl_test1_), secondary_domain(secdom) {
3361 if (brickname.size() == 0)
3362 brickname =
"Generic source term assembly brick";
3364 set_flags(brickname,
true ,
3368 directvarname = directvarname_; directdataname = directdataname_;
3373 static bool check_compatibility_vl_test(model &md,
3374 const model::varnamelist vl_test) {
3375 model::varnamelist org;
3376 for (
size_type i = 0; i < vl_test.size(); ++i) {
3377 if (md.is_affine_dependent_variable(vl_test[i]))
3378 org.push_back(md.org_variable(vl_test[i]));
3380 for (
size_type i = 0; i < vl_test.size(); ++i)
3381 for (
size_type j = 0; j < org.size(); ++j)
3382 if (vl_test[i].compare(org[j]) == 0)
return false;
3387 (model &md,
const mesh_im &mim,
const std::string &expr,
size_type region,
3388 const std::string &brickname, std::string directvarname,
3389 const std::string &directdataname,
bool return_if_nonlin,
3390 const std::string &secondary_domain) {
3392 ga_workspace workspace(md);
3393 size_type order = workspace.add_expression(expr, mim, region, 1,
3395 GMM_ASSERT1(order <= 1,
"Wrong order for a source term");
3396 model::varnamelist vl, vl_test1, vl_test2, dl;
3397 bool is_lin = workspace.used_variables(vl, vl_test1, vl_test2, dl, 1);
3398 if (!is_lin && return_if_nonlin)
return size_type(-1);
3399 GMM_ASSERT1(is_lin,
"Nonlinear term");
3400 GMM_ASSERT1(check_compatibility_vl_test(md, vl_test1),
3401 "This brick do not support the assembly on both an affine "
3402 "dependent variable and its original variable. "
3403 "Split the brick.");
3405 if (directdataname.size()) {
3406 vl.push_back(directvarname);
3407 dl.push_back(directdataname);
3408 }
else directvarname =
"";
3410 pbrick pbr = std::make_shared<gen_source_term_assembly_brick>
3411 (expr, brickname, vl_test1, directvarname, directdataname,
3415 for (
size_type i = 0; i < vl_test1.size(); ++i)
3416 tl.push_back(model::term_description(vl_test1[i]));
3417 if (directdataname.size())
3418 tl.push_back(model::term_description(directvarname));
3420 return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
3425 const std::string &brickname,
const std::string &directvarname,
3426 const std::string &directdataname,
bool return_if_nonlin) {
3427 return add_source_term_(md, mim, expr, region, brickname, directvarname,
3428 directdataname, return_if_nonlin,
"");
3433 const std::string &secondary_domain,
3434 const std::string &brickname,
const std::string &directvarname,
3435 const std::string &directdataname,
bool return_if_nonlin) {
3436 return add_source_term_(md, mim, expr, region, brickname, directvarname,
3437 directdataname, return_if_nonlin, secondary_domain);
3446 struct gen_linear_assembly_brick :
public virtual_brick {
3450 model::varnamelist vl_test1, vl_test2;
3451 std::string secondary_domain;
3453 virtual void asm_real_tangent_terms(
const model &md,
size_type ib,
3454 const model::varnamelist &,
3455 const model::varnamelist &dl,
3456 const model::mimlist &mims,
3457 model::real_matlist &matl,
3458 model::real_veclist &,
3459 model::real_veclist &,
3461 build_version version)
const {
3462 GMM_ASSERT1(matl.size() == vl_test1.size(),
3463 "Wrong number of terms for Generic linear assembly brick");
3464 GMM_ASSERT1(mims.size() == 1,
3465 "Generic linear assembly brick needs one and only one "
3467 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0);
3468 for (
size_type i = 0; i < dl.size(); ++i)
3469 recompute_matrix = recompute_matrix ||
3470 md.is_var_newer_than_brick(dl[i], ib);
3472 if (recompute_matrix) {
3474 ga_workspace workspace(md, ga_workspace::inherit::ALL);
3475 workspace.add_expression(expr, *(mims[0]), region, 2, secondary_domain);
3476 GMM_TRACE2(name <<
": generic matrix assembly");
3477 workspace.assembly(2);
3478 const auto &R=workspace.assembled_matrix();
3479 for (
size_type i = 0; i < vl_test1.size(); ++i) {
3480 scalar_type alpha = scalar_type(1)
3481 / ( workspace.factor_of_variable(vl_test1[i]) *
3482 workspace.factor_of_variable(vl_test2[i]));
3483 const auto &I1=workspace.interval_of_variable(vl_test1[i]),
3484 &I2=workspace.interval_of_variable(vl_test2[i]);
3485 gmm::copy(gmm::scaled(gmm::sub_matrix(R, I1, I2), alpha),
3491 virtual std::string declare_volume_assembly_string
3492 (
const model &,
size_type,
const model::varnamelist &,
3493 const model::varnamelist &)
const {
3494 return is_lower_dim ? std::string() : expr;
3497 gen_linear_assembly_brick(
const std::string &expr_,
const mesh_im &mim,
3499 bool is_coer, std::string brickname,
3500 const model::varnamelist &vl_test1_,
3501 const model::varnamelist &vl_test2_,
3502 const std::string &secdom)
3503 : vl_test1(vl_test1_), vl_test2(vl_test2_), secondary_domain(secdom) {
3504 if (brickname.size() == 0) brickname =
"Generic linear assembly brick";
3506 is_lower_dim = mim.is_lower_dimensional();
3507 set_flags(brickname,
true ,
3514 static bool check_compatibility_vl_test(model &md,
3515 const model::varnamelist vl_test1,
3516 const model::varnamelist vl_test2) {
3517 for (
size_type i = 0; i < vl_test1.size(); ++i)
3518 for (
size_type j = 0; j < vl_test1.size(); ++j) {
3519 bool is1 = md.is_affine_dependent_variable(vl_test1[i]);
3520 bool is2 = md.is_affine_dependent_variable(vl_test2[i]);
3522 const std::string &org1
3523 = is1 ? md.org_variable(vl_test1[i]) : vl_test1[i];
3524 const std::string &org2
3525 = is2 ? md.org_variable(vl_test2[i]) : vl_test2[i];
3526 bool is1_bis = md.is_affine_dependent_variable(vl_test1[j]);
3527 bool is2_bis = md.is_affine_dependent_variable(vl_test2[j]);
3528 const std::string &org1_bis = is1_bis ? md.org_variable(vl_test1[j])
3530 const std::string &org2_bis = is2_bis ? md.org_variable(vl_test2[j])
3532 if (org1.compare(org1_bis) == 0 && org2.compare(org2_bis))
3542 (model &md,
const mesh_im &mim,
const std::string &expr,
size_type region,
3543 bool is_sym,
bool is_coercive,
const std::string &brickname,
3544 bool return_if_nonlin,
const std::string &secondary_domain) {
3546 ga_workspace workspace(md, ga_workspace::inherit::ALL);
3547 size_type order = workspace.add_expression(expr, mim, region,
3548 2, secondary_domain);
3549 model::varnamelist vl, vl_test1, vl_test2, dl;
3550 bool is_lin = workspace.used_variables(vl, vl_test1, vl_test2, dl, 2);
3552 if (!is_lin && return_if_nonlin)
return size_type(-1);
3553 GMM_ASSERT1(is_lin,
"Nonlinear term");
3554 if (order == 0) { is_coercive = is_sym =
true; }
3556 std::string const_expr= workspace.extract_constant_term(mim.linked_mesh());
3557 if (const_expr.size()) {
3559 (md, mim, const_expr, region, brickname+
" (source term)",
3560 "",
"",
false, secondary_domain);
3565 GMM_ASSERT1(check_compatibility_vl_test(md, vl_test1, vl_test2),
3566 "This brick do not support the assembly on both an affine "
3567 "dependent variable and its original variable. "
3568 "Split the brick.");
3570 if (vl_test1.size()) {
3571 pbrick pbr = std::make_shared<gen_linear_assembly_brick>
3572 (expr, mim, is_sym, is_coercive, brickname, vl_test1, vl_test2,
3575 for (
size_type i = 0; i < vl_test1.size(); ++i)
3576 tl.push_back(model::term_description(vl_test1[i], vl_test2[i],
false));
3578 return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
3585 bool is_sym,
bool is_coercive,
const std::string &brickname,
3586 bool return_if_nonlin) {
3587 return add_linear_term_(md, mim, expr, region, is_sym, is_coercive,
3588 brickname, return_if_nonlin,
"");
3593 const std::string &secondary_domain,
bool is_sym,
bool is_coercive,
3594 const std::string &brickname,
bool return_if_nonlin) {
3595 return add_linear_term_(md, mim, expr, region, is_sym, is_coercive,
3596 brickname, return_if_nonlin, secondary_domain);
3606 struct gen_nonlinear_assembly_brick :
public virtual_brick {
3610 std::string secondary_domain;
3612 virtual void real_post_assembly_in_serial(
const model &md,
size_type ,
3613 const model::varnamelist &,
3614 const model::varnamelist &,
3615 const model::mimlist &mims,
3616 model::real_matlist &,
3617 model::real_veclist &,
3618 model::real_veclist &,
3620 build_version)
const {
3621 GMM_ASSERT1(mims.size() == 1,
3622 "Generic linear assembly brick needs one and only one "
3624 md.add_generic_expression(expr, *(mims[0]), region, secondary_domain);
3627 virtual std::string declare_volume_assembly_string
3628 (
const model &,
size_type,
const model::varnamelist &,
3629 const model::varnamelist &)
const {
3634 gen_nonlinear_assembly_brick(
const std::string &expr_,
const mesh_im &mim,
3637 std::string brickname,
3638 const std::string &secdom) {
3639 if (brickname.size() == 0) brickname =
"Generic linear assembly brick";
3641 secondary_domain = secdom;
3642 is_lower_dim = mim.is_lower_dimensional();
3643 set_flags(brickname,
false ,
3651 (model &md,
const mesh_im &mim,
const std::string &expr,
size_type region,
3652 bool is_sym,
bool is_coercive,
const std::string &brickname,
3653 const std::string &secondary_domain) {
3655 ga_workspace workspace(md);
3656 size_type order = workspace.add_expression(expr, mim, region, 2,
3658 GMM_ASSERT1(order < 2,
"Order two test functions (Test2) are not allowed"
3659 " in assembly string for nonlinear terms");
3660 model::varnamelist vl, vl_test1, vl_test2, ddl, dl;
3661 workspace.used_variables(vl, vl_test1, vl_test2, ddl, order);
3662 for (
size_type i = 0; i < ddl.size(); ++i)
3663 if (md.is_true_data(ddl[i])) dl.push_back(ddl[i]);
3664 else vl.push_back(ddl[i]);
3665 if (order == 0) { is_coercive = is_sym =
true; }
3666 pbrick pbr = std::make_shared<gen_nonlinear_assembly_brick>
3667 (expr, mim, is_sym, is_coercive, brickname, secondary_domain);
3671 return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
3677 bool is_sym,
bool is_coercive,
const std::string &brickname) {
3678 return add_nonlinear_term_(md, mim, expr, region, is_sym, is_coercive,
3684 const std::string &secondary_domain,
bool is_sym,
bool is_coercive,
3685 const std::string &brickname) {
3686 return add_nonlinear_term_(md, mim, expr, region, is_sym, is_coercive,
3687 brickname, secondary_domain);
3698 struct generic_elliptic_brick :
public virtual_brick {
3700 virtual void asm_real_tangent_terms(
const model &md,
size_type ,
3701 const model::varnamelist &vl,
3702 const model::varnamelist &dl,
3703 const model::mimlist &mims,
3704 model::real_matlist &matl,
3705 model::real_veclist &,
3706 model::real_veclist &,
3708 build_version)
const {
3709 GMM_ASSERT1(matl.size() == 1,
3710 "Generic elliptic brick has one and only one term");
3711 GMM_ASSERT1(mims.size() == 1,
3712 "Generic elliptic brick need one and only one mesh_im");
3713 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
3714 "Wrong number of variables for generic elliptic brick");
3716 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
3717 const mesh &m = mf_u.linked_mesh();
3718 size_type N = m.dim(), Q = mf_u.get_qdim(), s = 1;
3719 const mesh_im &mim = *mims[0];
3720 const model_real_plain_vector *A = 0;
3721 const mesh_fem *mf_a = 0;
3722 mesh_region rg(region);
3723 m.intersect_with_mpi_region(rg);
3724 if (dl.size() > 0) {
3725 A = &(md.real_variable(dl[0]));
3726 mf_a = md.pmesh_fem_of_variable(dl[0]);
3727 s = gmm::vect_size(*A);
3728 if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
3732 GMM_TRACE2(
"Generic elliptic term assembly");
3737 (matl[0], mim, mf_u, *mf_a, *A, rg);
3740 (matl[0], mim, mf_u, *mf_a, *A, rg);
3745 (matl[0], mim, mf_u, rg);
3748 (matl[0], mim, mf_u, rg);
3749 if (A) gmm::scale(matl[0], (*A)[0]);
3751 }
else if (s == N*N) {
3755 (matl[0], mim, mf_u, *mf_a, *A, rg);
3758 (matl[0], mim, mf_u, *mf_a, *A, rg);
3762 (matl[0], mim, mf_u, *A, rg);
3765 (matl[0], mim, mf_u, *A, rg);
3767 }
else if (s == N*N*Q*Q) {
3770 (matl[0], mim, mf_u, *mf_a, *A, rg);
3773 (matl[0], mim, mf_u, *A, rg);
3775 GMM_ASSERT1(
false,
"Bad format generic elliptic brick coefficient");
3779 virtual void real_post_assembly_in_serial(
const model &md,
size_type,
3780 const model::varnamelist &,
3781 const model::varnamelist &dl,
3782 const model::mimlist &,
3783 model::real_matlist &,
3784 model::real_veclist &,
3785 model::real_veclist &,
3787 build_version)
const
3789 const model_real_plain_vector *A = 0;
3790 const mesh_fem *mf_a = 0;
3793 A = &(md.real_variable(dl[0]));
3794 mf_a = md.pmesh_fem_of_variable(dl[0]);
3798 virtual void complex_post_assembly_in_serial(
const model &md,
size_type,
3799 const model::varnamelist &,
3800 const model::varnamelist &dl,
3801 const model::mimlist &,
3802 model::complex_matlist &,
3803 model::complex_veclist &,
3804 model::complex_veclist &,
3806 build_version)
const
3808 const model_real_plain_vector *A = 0;
3809 const mesh_fem *mf_a = 0;
3812 A = &(md.real_variable(dl[0]));
3813 mf_a = md.pmesh_fem_of_variable(dl[0]);
3817 virtual scalar_type asm_complex_tangent_terms(
const model &md,
size_type,
3818 const model::varnamelist &vl,
3819 const model::varnamelist &,
3820 const model::mimlist &,
3821 model::complex_matlist &matl,
3822 model::complex_veclist &,
3823 model::complex_veclist &,
3825 const model_complex_plain_vector &U = md.complex_variable(vl[0]);
3826 return gmm::abs(gmm::vect_hp(matl[0], U, U)) / scalar_type(2);
3830 virtual void asm_complex_tangent_terms(
const model &md,
size_type,
3831 const model::varnamelist &vl,
3832 const model::varnamelist &dl,
3833 const model::mimlist &mims,
3834 model::complex_matlist &matl,
3835 model::complex_veclist &,
3836 model::complex_veclist &,
3838 build_version)
const {
3839 GMM_ASSERT1(matl.size() == 1,
3840 "Generic elliptic brick has one and only one term");
3841 GMM_ASSERT1(mims.size() == 1,
3842 "Generic elliptic brick need one and only one mesh_im");
3843 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
3844 "Wrong number of variables for generic elliptic brick");
3846 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
3847 const mesh &m = mf_u.linked_mesh();
3848 size_type N = m.dim(), Q = mf_u.get_qdim(), s = 1;
3849 const mesh_im &mim = *mims[0];
3850 const model_real_plain_vector *A = 0;
3851 const mesh_fem *mf_a = 0;
3852 mesh_region rg(region);
3853 m.intersect_with_mpi_region(rg);
3856 if (dl.size() > 0) {
3857 A = &(md.real_variable(dl[0]));
3858 mf_a = md.pmesh_fem_of_variable(dl[0]);
3859 s = gmm::vect_size(*A);
3860 if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
3864 GMM_TRACE2(
"Generic elliptic term assembly");
3869 (matl[0], mim, mf_u, *mf_a, *A, rg);
3872 (matl[0], mim, mf_u, *mf_a, *A, rg);
3877 (gmm::real_part(matl[0]), mim, mf_u, rg);
3880 (gmm::real_part(matl[0]), mim, mf_u, rg);
3881 if (A) gmm::scale(matl[0], (*A)[0]);
3883 }
else if (s == N*N) {
3887 (matl[0], mim, mf_u, *mf_a, *A, rg);
3890 (matl[0], mim, mf_u, *mf_a, *A, rg);
3894 (matl[0], mim, mf_u, *A, rg);
3897 (matl[0], mim, mf_u, *A, rg);
3899 }
else if (s == N*N*Q*Q) {
3902 (matl[0], mim, mf_u, *mf_a, *A, rg);
3905 (matl[0], mim, mf_u, *A, rg);
3908 "Bad format generic elliptic brick coefficient");
3911 generic_elliptic_brick() {
3912 set_flags(
"Generic elliptic",
true ,
3920 const std::string &varname,
3923 pbrick pbr = std::make_shared<generic_elliptic_brick>();
3925 tl.push_back(model::term_description(varname, varname,
true));
3926 return md.
add_brick(pbr, model::varnamelist(1, varname),
3927 model::varnamelist(), tl, model::mimlist(1, &mim),
3930 std::string test_varname
3931 =
"Test_" + sup_previous_and_dot_to_varname(varname);
3936 expr =
"Grad_"+varname+
".Grad_"+test_varname;
3938 expr =
"Grad_"+varname+
":Grad_"+test_varname;
3940 "Laplacian",
false);
3945 const std::string &varname,
3946 const std::string &dataname,
3949 pbrick pbr = std::make_shared<generic_elliptic_brick>();
3951 tl.push_back(model::term_description(varname, varname,
true));
3952 return md.
add_brick(pbr, model::varnamelist(1, varname),
3953 model::varnamelist(1, dataname), tl,
3954 model::mimlist(1, &mim), region);
3956 std::string test_varname
3957 =
"Test_" + sup_previous_and_dot_to_varname(varname);
3970 if (qdim_data != 1) {
3971 GMM_ASSERT1(qdim_data == gmm::sqr(dim),
3972 "Wrong data size for generic elliptic brick");
3973 expr =
"((Reshape("+dataname+
",meshdim,meshdim))*Grad_"+varname+
").Grad_"
3976 expr =
"(("+dataname+
")*Grad_"+varname+
").Grad_"+test_varname;
3979 if (qdim_data != 1) {
3980 if (qdim_data == gmm::sqr(dim))
3981 expr =
"((Reshape("+dataname+
",meshdim,meshdim))*Grad_"+varname+
"):Grad_"
3983 else if (qdim_data == gmm::sqr(gmm::sqr(dim)))
3984 expr =
"((Reshape("+dataname+
",meshdim,meshdim,meshdim,meshdim))*Grad_"
3985 +varname+
"):Grad_"+test_varname;
3986 else GMM_ASSERT1(
false,
"Wrong data size for generic elliptic brick");
3988 expr =
"(("+dataname+
")*Grad_"+varname+
"):Grad_"+test_varname;
3992 (md, mim, expr, region,
true,
true,
"Generic elliptic",
true);
3995 "Generic elliptic (nonlinear)");
4007 struct source_term_brick :
public virtual_brick {
4009 void asm_real_tangent_terms(
const model &md,
size_type ,
4010 const model::varnamelist &vl,
4011 const model::varnamelist &dl,
4012 const model::mimlist &mims,
4013 model::real_matlist &,
4014 model::real_veclist &vecl,
4015 model::real_veclist &,
4017 build_version)
const override {
4018 GMM_ASSERT1(vecl.size() == 1,
4019 "Source term brick has one and only one term");
4020 GMM_ASSERT1(mims.size() == 1,
4021 "Source term brick need one and only one mesh_im");
4022 GMM_ASSERT1(vl.size() == 1 && dl.size() > 0 && dl.size() <= 2,
4023 "Wrong number of variables for source term brick");
4025 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4026 const mesh_im &mim = *mims[0];
4027 const model_real_plain_vector &A = md.real_variable(dl[0]);
4028 const mesh_fem *mf_data = md.pmesh_fem_of_variable(dl[0]);
4029 mesh_region rg(region);
4030 mim.linked_mesh().intersect_with_mpi_region(rg);
4033 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4035 GMM_ASSERT1(mf_u.get_qdim() == s,
4036 dl[0] <<
": bad format of source term data. "
4037 "Detected dimension is " << s <<
" should be "
4040 GMM_TRACE2(
"Source term assembly");
4046 if (dl.size() > 1) gmm::add(md.real_variable(dl[1]), vecl[0]);
4049 void real_post_assembly_in_serial(
const model &md,
size_type ib,
4050 const model::varnamelist &,
4051 const model::varnamelist &,
4052 const model::mimlist &,
4053 model::real_matlist &,
4054 model::real_veclist &vecl,
4055 model::real_veclist &,
4057 build_version)
const override
4059 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4063 void asm_complex_tangent_terms(
const model &md,
size_type ,
4064 const model::varnamelist &vl,
4065 const model::varnamelist &dl,
4066 const model::mimlist &mims,
4067 model::complex_matlist &,
4068 model::complex_veclist &vecl,
4069 model::complex_veclist &,
4071 build_version)
const override {
4072 GMM_ASSERT1(vecl.size() == 1,
4073 "Source term brick has one and only one term");
4074 GMM_ASSERT1(mims.size() == 1,
4075 "Source term brick need one and only one mesh_im");
4076 GMM_ASSERT1(vl.size() == 1 && dl.size() > 0 && dl.size() <= 2,
4077 "Wrong number of variables for source term brick");
4079 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4080 const mesh_im &mim = *mims[0];
4081 const model_complex_plain_vector &A = md.complex_variable(dl[0]);
4082 const mesh_fem *mf_data = md.pmesh_fem_of_variable(dl[0]);
4083 mesh_region rg(region);
4084 mim.linked_mesh().intersect_with_mpi_region(rg);
4087 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4089 GMM_ASSERT1(mf_u.get_qdim() == s,
"Bad format of source term data");
4091 GMM_TRACE2(
"Source term assembly");
4097 if (dl.size() > 1)
gmm::add(md.complex_variable(dl[1]), vecl[0]);
4100 void complex_post_assembly_in_serial(
const model &md,
4102 const model::varnamelist &,
4103 const model::varnamelist &,
4104 const model::mimlist &,
4105 model::complex_matlist &,
4106 model::complex_veclist &vecl,
4107 model::complex_veclist &,
4108 size_type, build_version)
const override
4110 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4115 source_term_brick() {
4116 set_flags(
"Source term",
true ,
4126 const std::string &varname,
4127 const std::string &dataexpr,
4129 const std::string &directdataname) {
4131 pbrick pbr = std::make_shared<source_term_brick>();
4133 tl.push_back(model::term_description(varname));
4134 model::varnamelist vdata(1, dataexpr);
4135 if (directdataname.size()) vdata.push_back(directdataname);
4136 return md.
add_brick(pbr, model::varnamelist(1, varname),
4137 vdata, tl, model::mimlist(1, &mim), region);
4139 std::string test_varname
4140 =
"Test_" + sup_previous_and_dot_to_varname(varname);
4145 expr =
"("+dataexpr+
")*"+test_varname;
4147 expr =
"("+dataexpr+
")."+test_varname;
4148 size_type ib = add_source_term_generic_assembly_brick
4149 (md, mim, expr, region,
"Source term", varname, directdataname,
true);
4152 "Source term (nonlinear)");
4153 if (directdataname.size())
4154 add_source_term_generic_assembly_brick
4155 (md, mim,
"", region,
"Source term", varname, directdataname);
4167 struct normal_source_term_brick :
public virtual_brick {
4169 void asm_real_tangent_terms(
const model &md,
size_type ,
4170 const model::varnamelist &vl,
4171 const model::varnamelist &dl,
4172 const model::mimlist &mims,
4173 model::real_matlist &,
4174 model::real_veclist &vecl,
4175 model::real_veclist &,
4177 build_version)
const override {
4178 GMM_ASSERT1(vecl.size() == 1,
4179 "Source term brick has one and only one term");
4180 GMM_ASSERT1(mims.size() == 1,
4181 "Source term brick need one and only one mesh_im");
4182 GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
4183 "Wrong number of variables for source term brick");
4185 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4186 const mesh_im &mim = *mims[0];
4187 const model_real_plain_vector &A = md.real_variable(dl[0]);
4188 const mesh_fem *mf_data = md.pmesh_fem_of_variable(dl[0]);
4189 mesh_region rg(region);
4190 mim.linked_mesh().intersect_with_mpi_region(rg);
4192 size_type s = gmm::vect_size(A), N = mf_u.linked_mesh().dim();
4193 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4195 GMM_ASSERT1(mf_u.get_qdim()*N == s,
4196 dl[0] <<
": bad format of normal source term data. "
4197 "Detected dimension is " << s <<
" should be "
4200 GMM_TRACE2(
"source term assembly");
4207 void real_post_assembly_in_serial(
const model &md,
size_type ib,
4208 const model::varnamelist &,
4209 const model::varnamelist &,
4210 const model::mimlist &,
4211 model::real_matlist &,
4212 model::real_veclist &vecl,
4213 model::real_veclist &,
4215 build_version)
const override {
4216 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4220 virtual void asm_complex_tangent_terms(
const model &md,
size_type ,
4221 const model::varnamelist &vl,
4222 const model::varnamelist &dl,
4223 const model::mimlist &mims,
4224 model::complex_matlist &,
4225 model::complex_veclist &vecl,
4226 model::complex_veclist &,
4228 build_version)
const {
4229 GMM_ASSERT1(vecl.size() == 1,
4230 "Source term brick has one and only one term");
4231 GMM_ASSERT1(mims.size() == 1,
4232 "Source term brick need one and only one mesh_im");
4233 GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
4234 "Wrong number of variables for source term brick");
4236 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4237 const mesh_im &mim = *mims[0];
4238 const model_complex_plain_vector &A = md.complex_variable(dl[0]);
4239 const mesh_fem *mf_data = md.pmesh_fem_of_variable(dl[0]);
4240 mesh_region rg(region);
4241 mim.linked_mesh().intersect_with_mpi_region(rg);
4243 size_type s = gmm::vect_size(A), N = mf_u.linked_mesh().dim();
4244 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4246 GMM_ASSERT1(s == mf_u.get_qdim()*N,
"Bad format of source term data");
4248 GMM_TRACE2(
"source term assembly");
4255 void complex_post_assembly_in_serial(
const model &md,
size_type ib,
4256 const model::varnamelist &,
4257 const model::varnamelist &,
4258 const model::mimlist &,
4259 model::complex_matlist &,
4260 model::complex_veclist &vecl,
4261 model::complex_veclist &,
4263 build_version)
const override {
4264 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4267 normal_source_term_brick() {
4268 set_flags(
"Normal source term",
true ,
4278 const std::string &varname,
4279 const std::string &dataexpr,
4282 pbrick pbr = std::make_shared<normal_source_term_brick>();
4284 tl.push_back(model::term_description(varname));
4285 model::varnamelist vdata(1, dataexpr);
4286 return md.
add_brick(pbr, model::varnamelist(1, varname),
4287 vdata, tl, model::mimlist(1, &mim), region);
4289 std::string test_varname
4290 =
"Test_" + sup_previous_and_dot_to_varname(varname);
4295 expr =
"(("+dataexpr+
").Normal)*"+test_varname;
4297 expr =
"(Reshape("+dataexpr+
",qdim("+varname
4298 +
"),meshdim)*Normal)."+test_varname;
4299 return add_source_term_generic_assembly_brick
4300 (md, mim, expr, region,
"Source term");
4313 struct Dirichlet_condition_brick :
public virtual_brick {
4316 bool normal_component;
4317 const mesh_fem *mf_mult_;
4323 virtual void asm_real_tangent_terms(
const model &md,
size_type ib,
4324 const model::varnamelist &vl,
4325 const model::varnamelist &dl,
4326 const model::mimlist &mims,
4327 model::real_matlist &matl,
4328 model::real_veclist &vecl,
4329 model::real_veclist &,
4331 build_version version)
const {
4332 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
4333 "Dirichlet condition brick has one and only one term");
4334 GMM_ASSERT1(mims.size() == 1,
4335 "Dirichlet condition brick need one and only one mesh_im");
4336 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() <= 3,
4337 "Wrong number of variables for Dirichlet condition brick");
4339 model_real_sparse_matrix& rB = rB_th;
4340 model_real_plain_vector& rV = rV_th;
4342 bool penalized = (vl.size() == 1);
4343 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4344 const mesh_fem &mf_mult = penalized ? (mf_mult_ ? *mf_mult_ : mf_u)
4345 : md.mesh_fem_of_variable(vl[1]);
4346 const mesh_im &mim = *mims[0];
4347 const model_real_plain_vector *A = 0, *COEFF = 0, *H = 0;
4348 const mesh_fem *mf_data = 0, *mf_H = 0;
4349 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
4350 || (penalized && md.is_var_newer_than_brick(dl[0], ib));
4353 COEFF = &(md.real_variable(dl[0]));
4354 GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
4355 "Data for coefficient should be a scalar");
4358 size_type s = 0, ind = (penalized ? 1 : 0);
4359 if (dl.size() > ind) {
4360 A = &(md.real_variable(dl[ind]));
4361 mf_data = md.pmesh_fem_of_variable(dl[ind]);
4362 s = gmm::vect_size(*A);
4363 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4364 size_type ss = ((normal_component) ? 1 : mf_u.get_qdim());
4365 GMM_ASSERT1(s == ss, dl[ind] <<
": bad format of "
4366 "Dirichlet data. Detected dimension is " << s
4367 <<
" should be " << ss);
4370 if (dl.size() > ind + 1) {
4371 GMM_ASSERT1(H_version,
4372 "Wrong number of data for Dirichlet condition brick");
4373 H = &(md.real_variable(dl[ind+1]));
4374 mf_H = md.pmesh_fem_of_variable(dl[ind+1]);
4375 s = gmm::vect_size(*A);
4377 s = s * mf_H->get_qdim() / mf_H->nb_dof();
4378 GMM_ASSERT1(mf_H->get_qdim() == 1,
"Implemented only for mf_H "
4379 "a scalar finite element method");
4381 GMM_ASSERT1(s = gmm::sqr(mf_u.get_qdim()),
4382 dl[ind+1] <<
": bad format of Dirichlet data. "
4383 "Detected dimension is " << s <<
" should be "
4384 <<
size_type(gmm::sqr(mf_u.get_qdim())));
4387 mesh_region rg(region);
4388 mim.linked_mesh().intersect_with_mpi_region(rg);
4390 if (recompute_matrix) {
4391 model_real_sparse_matrix *B = &(matl[0]);
4392 if (penalized && (&mf_mult != &mf_u)) {
4399 GMM_TRACE2(
"Mass term assembly for Dirichlet condition");
4400 if (H_version || normal_component) {
4401 ga_workspace workspace;
4402 gmm::sub_interval Imult(0, mf_mult.nb_dof()), Iu(0, mf_u.nb_dof());
4403 base_vector u(mf_u.nb_dof());
4404 base_vector
mult(mf_mult.nb_dof());
4405 workspace.add_fem_variable(
"u", mf_u, Iu, u);
4406 workspace.add_fem_variable(
"mult", mf_mult, Imult, mult);
4407 auto expression = std::string{};
4409 if (mf_H) workspace.add_fem_constant(
"A", *mf_H, *H);
4410 else workspace.add_fixed_size_constant(
"A", *H);
4411 expression = (mf_u.get_qdim() == 1) ?
4412 "Test_mult . (A . Test2_u)"
4414 "Test_mult. (Reshape(A, qdim(u), qdim(u)) . Test2_u)";
4415 }
else if (normal_component) {
4416 expression =
"Test_mult . (Test2_u . Normal)";
4418 workspace.add_expression(expression, mim, rg);
4419 workspace.set_assembled_matrix(*B);
4420 workspace.assembly(2);
4425 if (penalized && (&mf_mult != &mf_u)) {
4426 GMM_ASSERT1(MPI_IS_MASTER(),
"Sorry, the penalized option "
4427 "filtered by a multiplier space is not parallelized");
4428 gmm::mult(gmm::transposed(rB), rB, matl[0]);
4429 gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
4430 }
else if (penalized) {
4431 gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
4435 if (dl.size() > ind) {
4436 GMM_TRACE2(
"Source term assembly for Dirichlet condition");
4438 if (penalized && (&mf_mult != &mf_u)) {
4452 if (penalized && (&mf_mult != &mf_u)) {
4453 gmm::mult(gmm::transposed(rB), rV, vecl[0]);
4454 gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
4455 rV = model_real_plain_vector();
4456 }
else if (penalized)
4457 gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
4462 void real_post_assembly_in_serial(
const model &md,
size_type ib,
4463 const model::varnamelist &,
4464 const model::varnamelist &,
4465 const model::mimlist &,
4466 model::real_matlist &,
4467 model::real_veclist &vecl,
4468 model::real_veclist &,
4470 build_version)
const override {
4471 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4474 virtual void asm_complex_tangent_terms(
const model &md,
size_type ib,
4475 const model::varnamelist &vl,
4476 const model::varnamelist &dl,
4477 const model::mimlist &mims,
4478 model::complex_matlist &matl,
4479 model::complex_veclist &vecl,
4480 model::complex_veclist &,
4482 build_version version)
const {
4483 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
4484 "Dirichlet condition brick has one and only one term");
4485 GMM_ASSERT1(mims.size() == 1,
4486 "Dirichlet condition brick need one and only one mesh_im");
4487 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() <= 3,
4488 "Wrong number of variables for Dirichlet condition brick");
4490 model_complex_sparse_matrix& cB = cB_th;
4491 model_complex_plain_vector& cV = cV_th;
4493 bool penalized = (vl.size() == 1);
4494 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4495 const mesh_fem &mf_mult = penalized ? (mf_mult_ ? *mf_mult_ : mf_u)
4496 : md.mesh_fem_of_variable(vl[1]);
4497 const mesh_im &mim = *mims[0];
4498 const model_complex_plain_vector *A = 0, *COEFF = 0, *H = 0;
4499 const mesh_fem *mf_data = 0, *mf_H = 0;
4500 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
4501 || (penalized && md.is_var_newer_than_brick(dl[0], ib));
4504 COEFF = &(md.complex_variable(dl[0]));
4505 GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
4506 "Data for coefficient should be a scalar");
4509 size_type s = 0, ind = (penalized ? 1 : 0);
4510 if (dl.size() > ind) {
4511 A = &(md.complex_variable(dl[ind]));
4512 mf_data = md.pmesh_fem_of_variable(dl[ind]);
4513 s = gmm::vect_size(*A);
4514 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4515 size_type ss = s * ((normal_component) ? mf_u.linked_mesh().dim() : 1);
4516 GMM_ASSERT1(mf_u.get_qdim() == ss,
4517 dl[ind] <<
": bad format of Dirichlet data. "
4518 "Detected dimension is " << ss <<
" should be "
4522 if (dl.size() > ind + 1) {
4523 GMM_ASSERT1(H_version,
4524 "Wrong number of data for Dirichlet condition brick");
4525 H = &(md.complex_variable(dl[ind+1]));
4526 mf_H = md.pmesh_fem_of_variable(dl[ind+1]);
4527 s = gmm::vect_size(*A);
4529 s = s * mf_H->get_qdim() / mf_H->nb_dof();
4530 GMM_ASSERT1(mf_H->get_qdim() == 1,
"Implemented only for mf_H "
4531 "a scalar finite element method");
4533 GMM_ASSERT1(s = gmm::sqr(mf_u.get_qdim()),
4534 dl[ind+1] <<
": bad format of Dirichlet data. "
4535 "Detected dimension is " << s <<
" should be "
4536 <<
size_type(gmm::sqr(mf_u.get_qdim())));
4539 mesh_region rg(region);
4540 mim.linked_mesh().intersect_with_mpi_region(rg);
4542 if (recompute_matrix) {
4543 model_complex_sparse_matrix *B = &(matl[0]);
4544 if (penalized && (&mf_mult != &mf_u)) {
4551 GMM_TRACE2(
"Mass term assembly for Dirichlet condition");
4553 if (mf_u.get_qdim() == 1)
4554 asm_real_or_complex_1_param_mat(*B, mim, mf_mult, mf_H, *H, rg,
4555 "(A*Test_u).Test2_u");
4557 asm_real_or_complex_1_param_mat(*B, mim, mf_mult, mf_H, *H, rg,
4558 "(Reshape(A,qdim(u),qdim(u))*Test2_u).Test_u");
4574 else if (normal_component) {
4575 ga_workspace workspace;
4576 gmm::sub_interval Imult(0, mf_mult.nb_dof()), Iu(0, mf_u.nb_dof());
4577 base_vector
mult(mf_mult.nb_dof()), u(mf_u.nb_dof());
4578 workspace.add_fem_variable(
"mult", mf_mult, Imult, mult);
4579 workspace.add_fem_variable(
"u", mf_u, Iu, u);
4580 workspace.add_expression(
"Test_mult.(Test2_u.Normal)", mim, rg);
4581 model_real_sparse_matrix BB(mf_mult.nb_dof(), mf_u.nb_dof());
4582 workspace.set_assembled_matrix(BB);
4583 workspace.assembly(2);
4599 if (penalized && (&mf_mult != &mf_u)) {
4600 gmm::mult(gmm::transposed(cB), cB, matl[0]);
4601 gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
4602 }
else if (penalized) {
4603 gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
4607 if (dl.size() > ind) {
4608 GMM_TRACE2(
"Source term assembly for Dirichlet condition");
4610 if (penalized && (&mf_mult != &mf_u)) {
4624 if (penalized && (&mf_mult != &mf_u)) {
4625 gmm::mult(gmm::transposed(cB), cV, vecl[0]);
4626 gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
4627 cV = model_complex_plain_vector();
4628 }
else if (penalized)
4629 gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
4633 void complex_post_assembly_in_serial(
const model &md,
size_type ib,
4634 const model::varnamelist &,
4635 const model::varnamelist &,
4636 const model::mimlist &,
4637 model::complex_matlist &,
4638 model::complex_veclist &vecl,
4639 model::complex_veclist &,
4641 build_version)
const override {
4642 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4646 virtual std::string declare_volume_assembly_string
4647 (
const model &,
size_type,
const model::varnamelist &,
4648 const model::varnamelist &)
const {
4649 return std::string();
4652 Dirichlet_condition_brick(
bool penalized,
bool H_version_,
4653 bool normal_component_,
4654 const mesh_fem *mf_mult__ = 0) {
4655 mf_mult_ = mf_mult__;
4656 H_version = H_version_;
4657 normal_component = normal_component_;
4658 GMM_ASSERT1(!(H_version && normal_component),
"Bad Dirichlet version");
4659 set_flags(penalized ?
"Dirichlet with penalization brick"
4660 :
"Dirichlet with multipliers brick",
4669 (
model &md,
const mesh_im &mim,
const std::string &varname,
4670 const std::string &multname,
size_type region,
4671 const std::string &dataname) {
4672 pbrick pbr = std::make_shared<Dirichlet_condition_brick>(
false,
false,
false);
4674 tl.push_back(model::term_description(multname, varname,
true));
4675 model::varnamelist vl(1, varname);
4676 vl.push_back(multname);
4677 model::varnamelist dl;
4678 if (dataname.size()) dl.push_back(dataname);
4679 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4683 (
model &md,
const mesh_im &mim,
const std::string &varname,
4685 const std::string &dataname) {
4686 std::string multname = md.
new_name(
"mult_on_" + varname);
4689 (md, mim, varname, multname, region, dataname);
4693 (
model &md,
const mesh_im &mim,
const std::string &varname,
4695 const std::string &dataname) {
4700 (md, mim, varname, mf_mult, region, dataname);
4708 (
model &md,
const mesh_im &mim,
const std::string &varname,
4709 scalar_type penalisation_coeff,
size_type region,
4710 const std::string &dataname,
const mesh_fem *mf_mult) {
4711 std::string coeffname = md.
new_name(
"penalization_on_" + varname);
4717 pbrick pbr = std::make_shared<Dirichlet_condition_brick>
4718 (
true,
false,
false, mf_mult);
4720 tl.push_back(model::term_description(varname, varname,
true));
4721 model::varnamelist vl(1, varname);
4722 model::varnamelist dl(1, coeffname);
4723 if (dataname.size()) dl.push_back(dataname);
4724 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4728 (
model &md,
const mesh_im &mim,
const std::string &varname,
4729 const std::string &multname,
size_type region,
4730 const std::string &dataname) {
4731 pbrick pbr = std::make_shared<Dirichlet_condition_brick>(
false,
false,
true);
4733 tl.push_back(model::term_description(multname, varname,
true));
4734 model::varnamelist vl(1, varname);
4735 vl.push_back(multname);
4736 model::varnamelist dl;
4737 if (dataname.size()) dl.push_back(dataname);
4738 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4742 (
model &md,
const mesh_im &mim,
const std::string &varname,
4744 const std::string &dataname) {
4745 std::string multname = md.
new_name(
"mult_on_" + varname);
4748 (md, mim, varname, multname, region, dataname);
4752 (
model &md,
const mesh_im &mim,
const std::string &varname,
4754 const std::string &dataname) {
4758 (md, mim, varname, mf_mult, region, dataname);
4762 (
model &md,
const mesh_im &mim,
const std::string &varname,
4763 scalar_type penalisation_coeff,
size_type region,
4764 const std::string &dataname,
const mesh_fem *mf_mult) {
4765 std::string coeffname = md.
new_name(
"penalization_on_" + varname);
4771 pbrick pbr = std::make_shared<Dirichlet_condition_brick>
4772 (
true,
false,
true, mf_mult);
4774 tl.push_back(model::term_description(varname, varname,
true));
4775 model::varnamelist vl(1, varname);
4776 model::varnamelist dl(1, coeffname);
4777 if (dataname.size()) dl.push_back(dataname);
4778 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4783 (
model &md,
const mesh_im &mim,
const std::string &varname,
4784 const std::string &multname,
size_type region,
4785 const std::string &dataname,
const std::string &Hname) {
4786 pbrick pbr = std::make_shared<Dirichlet_condition_brick>(
false,
true,
false);
4788 tl.push_back(model::term_description(multname, varname,
true));
4789 model::varnamelist vl(1, varname);
4790 vl.push_back(multname);
4791 model::varnamelist dl;
4792 dl.push_back(dataname);
4793 dl.push_back(Hname);
4794 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4798 (
model &md,
const mesh_im &mim,
const std::string &varname,
4800 const std::string &dataname,
const std::string &Hname) {
4801 std::string multname = md.
new_name(
"mult_on_" + varname);
4804 (md, mim, varname, multname, region, dataname, Hname);
4808 (
model &md,
const mesh_im &mim,
const std::string &varname,
4810 const std::string &dataname,
const std::string &Hname) {
4815 (md, mim, varname, mf_mult, region, dataname, Hname);
4819 (
model &md,
const mesh_im &mim,
const std::string &varname,
4820 scalar_type penalisation_coeff,
size_type region,
4821 const std::string &dataname,
const std::string &Hname,
4823 std::string coeffname = md.
new_name(
"penalization_on_" + varname);
4829 pbrick pbr = std::make_shared<Dirichlet_condition_brick>
4830 (
true,
true,
false, mf_mult);
4832 tl.push_back(model::term_description(varname, varname,
true));
4833 model::varnamelist vl(1, varname);
4834 model::varnamelist dl(1, coeffname);
4835 dl.push_back(dataname);
4836 dl.push_back(Hname);
4837 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4841 scalar_type penalisation_coeff) {
4845 GMM_ASSERT1(gmm::vect_size(d)==1,
4846 "Wrong coefficient size, may be not a Dirichlet brick "
4847 "with penalization");
4848 d[0] = penalisation_coeff;
4852 GMM_ASSERT1(gmm::vect_size(d)==1,
4853 "Wrong coefficient size, may be not a Dirichlet brick "
4854 "with penalization");
4855 d[0] = penalisation_coeff;
4865 struct simplification_Dirichlet_condition_brick :
public virtual_brick {
4867 virtual void real_pre_assembly_in_serial(
const model &md,
size_type ,
4868 const model::varnamelist &vl,
4869 const model::varnamelist &dl,
4870 const model::mimlist &mims,
4871 model::real_matlist &matl,
4872 model::real_veclist &vecl,
4873 model::real_veclist &,
4875 build_version )
const {
4876 if (MPI_IS_MASTER()) {
4878 GMM_ASSERT1(vecl.size() == 0 && matl.size() == 0,
4879 "Dirichlet condition brick by simplification has no term");
4880 GMM_ASSERT1(mims.size() == 0,
4881 "Dirichlet condition brick by simplification need no "
4883 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
4884 "Wrong number of variables for Dirichlet condition brick "
4885 "by simplification");
4887 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4888 const model_real_plain_vector *A = 0;
4889 const mesh_fem *mf_data = 0;
4892 if (dl.size() == 1) {
4893 A = &(md.real_variable(dl[0]));
4894 mf_data = md.pmesh_fem_of_variable(dl[0]);
4897 GMM_ASSERT1(mf_data == &mf_u,
"Sorry, for this brick, the data has"
4898 " to be defined on the same f.e.m. as the unknown");
4900 s = gmm::vect_size(*A);
4901 GMM_ASSERT1(mf_u.get_qdim() == s,
": bad format of "
4902 "Dirichlet data. Detected dimension is " << s
4903 <<
" should be " <<
size_type(mf_u.get_qdim()));
4907 mesh_region rg(region);
4910 if (mf_u.get_qdim() > 1 || (!mf_data && A)) {
4911 for (mr_visitor i(rg, mf_u.linked_mesh()); !i.finished(); ++i) {
4912 pfem pf = mf_u.fem_of_element(i.cv());
4914 GMM_ASSERT1(pf->target_dim() == 1,
4915 "Intrinsically vectorial fems are not allowed");
4916 GMM_ASSERT1(mf_data || pf->is_lagrange(),
"Constant Dirichlet "
4917 "data allowed for lagrange fems only");
4922 dal::bit_vector dofs = mf_u.dof_on_region(rg);
4924 if (A && !mf_data) {
4925 GMM_ASSERT1(dofs.card() % s == 0,
"Problem with dof vectorization");
4928 for (dal::bv_visitor i(dofs); !i.finished(); ++i) {
4930 if (A) val = (mf_data ? (*A)[i] : (*A)[i%s]);
4931 md.add_real_dof_constraint(vl[0], i, val);
4936 virtual void complex_pre_assembly_in_serial(
const model &md,
size_type ,
4937 const model::varnamelist &vl,
4938 const model::varnamelist &dl,
4939 const model::mimlist &mims,
4940 model::complex_matlist &matl,
4941 model::complex_veclist &vecl,
4942 model::complex_veclist &,
4944 build_version )
const {
4945 if (MPI_IS_MASTER()) {
4946 GMM_ASSERT1(vecl.size() == 0 && matl.size() == 0,
4947 "Dirichlet condition brick by simplification has no term");
4948 GMM_ASSERT1(mims.size() == 0,
4949 "Dirichlet condition brick by simplification need no "
4951 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
4952 "Wrong number of variables for Dirichlet condition brick "
4953 "by simplification");
4955 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4956 const model_complex_plain_vector *A = 0;
4957 const mesh_fem *mf_data = 0;
4960 if (dl.size() == 1) {
4961 A = &(md.complex_variable(dl[0]));
4962 mf_data = md.pmesh_fem_of_variable(dl[0]);
4965 GMM_ASSERT1(mf_data == &mf_u,
"Sorry, for this brick, the data has"
4966 " to be define on the same f.e.m. than the unknown");
4968 s = gmm::vect_size(*A);
4969 GMM_ASSERT1(mf_u.get_qdim() == s,
": bad format of "
4970 "Dirichlet data. Detected dimension is " << s
4971 <<
" should be " <<
size_type(mf_u.get_qdim()));
4975 mesh_region rg(region);
4978 if (mf_u.get_qdim() > 1 || (!mf_data && A)) {
4979 for (mr_visitor i(rg, mf_u.linked_mesh()); !i.finished(); ++i) {
4980 pfem pf = mf_u.fem_of_element(i.cv());
4982 GMM_ASSERT1(pf->target_dim() == 1,
4983 "Intrinsically vectorial fems are not allowed");
4984 GMM_ASSERT1(mf_data || pf->is_lagrange(),
"Constant Dirichlet "
4985 "data allowed for lagrange fems only");
4990 dal::bit_vector dofs = mf_u.dof_on_region(rg);
4992 if (A && !mf_data) {
4993 GMM_ASSERT1(dofs.card() % s == 0,
"Problem with dof vectorization");
4996 for (dal::bv_visitor i(dofs); !i.finished(); ++i) {
4997 complex_type val(0);
4998 if (A) val = (mf_data ? (*A)[i] : (*A)[i%s]);
4999 md.add_complex_dof_constraint(vl[0], i, val);
5005 virtual std::string declare_volume_assembly_string
5006 (
const model &,
size_type,
const model::varnamelist &,
5007 const model::varnamelist &)
const {
5008 return std::string();
5011 simplification_Dirichlet_condition_brick() {
5012 set_flags(
"Dirichlet with simplification brick",
5021 (
model &md,
const std::string &varname,
5022 size_type region,
const std::string &dataname) {
5023 pbrick pbr = std::make_shared<simplification_Dirichlet_condition_brick>();
5025 model::varnamelist vl(1, varname);
5026 model::varnamelist dl;
5027 if (dataname.size()) dl.push_back(dataname);
5028 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(), region);
5038 (
model &md,
const mesh_im &mim,
const std::string &varname,
5039 const std::string &Neumannterm,
5040 const std::string &datagamma0,
size_type region, scalar_type theta_,
5041 const std::string &datag) {
5042 std::string theta = std::to_string(theta_);
5044 ga_workspace workspace(md, ga_workspace::inherit::ALL);
5045 size_type order = workspace.add_expression(Neumannterm, mim, region, 1);
5046 GMM_ASSERT1(order == 0,
"Wrong expression of the Neumann term");
5047 bool is_lin = workspace.is_linear(1);
5049 std::string condition =
"("+varname + (datag.size() ?
"-("+datag+
"))":
")");
5050 std::string gamma =
"(("+datagamma0+
")*element_size)";
5051 std::string r =
"(1/"+gamma+
")";
5052 std::string expr =
"("+r+
"*"+condition+
"-("+Neumannterm+
")).Test_"+varname;
5053 if (theta_ != scalar_type(0)) {
5054 std::string derivative_Neumann = workspace.extract_order1_term(varname);
5055 if (derivative_Neumann.size())
5056 expr+=
"-"+theta+
"*"+condition+
".("+derivative_Neumann+
")";
5064 "Dirichlet condition with Nitsche's method");
5067 "Dirichlet condition with Nitsche's method");
5072 (
model &md,
const mesh_im &mim,
const std::string &varname,
5073 const std::string &Neumannterm,
5074 const std::string &datagamma0,
size_type region, scalar_type theta_,
5075 const std::string &datag) {
5076 std::string theta = std::to_string(theta_);
5078 ga_workspace workspace(md, ga_workspace::inherit::ALL);
5079 size_type order = workspace.add_expression(Neumannterm, mim, region, 1);
5080 GMM_ASSERT1(order == 0,
"Wrong expression of the Neumann term");
5081 bool is_lin = workspace.is_linear(1);
5083 std::string condition =
"("+varname+
".Normal"
5084 + (datag.size() ?
"-("+datag+
"))":
")");
5085 std::string gamma =
"(("+datagamma0+
")*element_size)";
5086 std::string r =
"(1/"+gamma+
")";
5087 std::string expr =
"("+r+
"*"+condition+
"-Normal.("+Neumannterm
5088 +
"))*(Normal.Test_"+varname+
")";
5089 if (theta_ != scalar_type(0)) {
5090 std::string derivative_Neumann = workspace.extract_order1_term(varname);
5091 if (derivative_Neumann.size())
5092 expr+=
"-"+theta+
"*"+condition+
"*Normal.("
5093 +derivative_Neumann+
")";
5097 "Dirichlet condition with Nitsche's method");
5100 "Dirichlet condition with Nitsche's method");
5105 (
model &md,
const mesh_im &mim,
const std::string &varname,
5106 const std::string &Neumannterm,
5107 const std::string &datagamma0,
size_type region, scalar_type theta_,
5108 const std::string &datag,
const std::string &dataH) {
5109 std::string theta = std::to_string(theta_);
5111 ga_workspace workspace(md, ga_workspace::inherit::ALL);
5112 size_type order = workspace.add_expression(Neumannterm, mim, region, 1);
5113 GMM_ASSERT1(order == 0,
"Wrong expression of the Neumann term");
5114 bool is_lin = workspace.is_linear(1);
5116 std::string condition =
"(("+dataH+
")*"+varname
5117 + (datag.size() ?
"-("+datag+
"))":
")");
5118 std::string gamma =
"(("+datagamma0+
")*element_size)";
5119 std::string r =
"(1/"+gamma+
")";
5120 std::string expr =
"("+r+
"*"+condition+
"-("+dataH+
")*("+Neumannterm
5121 +
"))*(("+dataH+
")*Test_"+varname+
")";
5122 if (theta_ != scalar_type(0)) {
5123 std::string derivative_Neumann = workspace.extract_order1_term(varname);
5124 if (derivative_Neumann.size())
5125 expr+=
"-"+theta+
"*"+condition+
"*(("+dataH+
")*("
5126 +derivative_Neumann+
"))";
5130 "Dirichlet condition with Nitsche's method");
5133 "Dirichlet condition with Nitsche's method");
5145 struct pointwise_constraints_brick :
public virtual_brick {
5147 mutable gmm::row_matrix<model_real_sparse_vector> rB;
5148 mutable gmm::row_matrix<model_complex_sparse_vector> cB;
5150 virtual void real_pre_assembly_in_serial(
const model &md,
size_type ib,
5151 const model::varnamelist &vl,
5152 const model::varnamelist &dl,
5153 const model::mimlist &mims,
5154 model::real_matlist &matl,
5155 model::real_veclist &vecl,
5156 model::real_veclist &,
5158 build_version version)
const {
5159 if (MPI_IS_MASTER()) {
5161 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5162 "Pointwize constraints brick has only one term");
5163 GMM_ASSERT1(mims.size() == 0,
5164 "Pointwize constraints brick does not need a mesh_im");
5165 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2,
5166 "Wrong number of variables for pointwize constraints brick");
5167 bool penalized = (vl.size() == 1);
5168 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5169 dim_type N = mf_u.linked_mesh().dim(), Q = mf_u.get_qdim(), ind_pt = 0;
5171 GMM_ASSERT1(dl.size() == dlsize || dl.size() == dlsize+1,
5172 "Wrong number of data for pointwize constraints brick");
5175 const model_real_plain_vector *COEFF = 0;
5177 COEFF = &(md.real_variable(dl[0]));
5179 GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
5180 "Data for coefficient should be a scalar");
5183 const model_real_plain_vector &PT = md.real_variable(dl[ind_pt]);
5184 size_type nb_co = gmm::vect_size(PT) / N;
5186 dim_type ind_unitv = dim_type((Q > 1) ? ind_pt+1 : 0);
5187 const model_real_plain_vector &unitv =md.real_variable(dl[ind_unitv]);
5188 GMM_ASSERT1((!ind_unitv || gmm::vect_size(unitv) == nb_co * Q),
5189 "Wrong size for vector of unit vectors");
5191 dim_type ind_rhs = dim_type((Q > 1) ? ind_pt+2 : ind_pt+1);
5192 if (dl.size() <
size_type(ind_rhs + 1)) ind_rhs = 0;
5193 const model_real_plain_vector &rhs = md.real_variable(dl[ind_rhs]);
5194 GMM_ASSERT1((!ind_rhs || gmm::vect_size(rhs) == nb_co),
5195 "Wrong size for vector of rhs");
5197 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
5198 || (penalized && (md.is_var_newer_than_brick(dl[ind_pt], ib)
5199 || md.is_var_newer_than_brick(dl[ind_unitv], ib)
5200 || md.is_var_newer_than_brick(dl[ind_rhs], ib)));
5202 if (recompute_matrix) {
5203 gmm::row_matrix<model_real_sparse_vector> BB(nb_co*Q, mf_u.nb_dof());
5206 dal::bit_vector dof_untouched;
5207 getfem::mesh_trans_inv mti(mf_u.linked_mesh());
5210 gmm::copy(gmm::sub_vector(PT, gmm::sub_interval(i*N, N)), pt);
5213 gmm::row_matrix<model_real_sparse_vector> &BBB = ((Q > 1) ? BB : rB);
5214 model_real_plain_vector vv;
5215 interpolation(mf_u, mti, vv, vv, BBB, 1, 1, &dof_untouched);
5216 GMM_ASSERT1(dof_untouched.card() == 0,
5217 "Pointwize constraints : some of the points are outside "
5218 "the mesh: " << dof_untouched);
5223 gmm::add(gmm::scaled(gmm::mat_row(BB, i*Q+q), unitv[i*Q+q]),
5224 gmm::mat_row(rB, i));
5227 gmm::mult(gmm::transposed(rB), rB, matl[0]);
5228 gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
5235 gmm::mult(gmm::transposed(rB), rhs, vecl[0]);
5236 gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
5244 virtual void complex_pre_assembly_in_serial(
const model &md,
size_type ib,
5245 const model::varnamelist &vl,
5246 const model::varnamelist &dl,
5247 const model::mimlist &mims,
5248 model::complex_matlist &matl,
5249 model::complex_veclist &vecl,
5250 model::complex_veclist &,
5252 build_version version)
const {
5253 if (MPI_IS_MASTER()) {
5254 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5255 "Pointwize constraints brick only one term");
5256 GMM_ASSERT1(mims.size() == 0,
5257 "Pointwize constraints brick does not need a mesh_im");
5258 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2,
5259 "Wrong number of variables for pointwize constraints brick");
5260 bool penalized = (vl.size() == 1);
5261 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5262 dim_type N = mf_u.linked_mesh().dim(), Q = mf_u.get_qdim(), ind_pt = 0;
5264 GMM_ASSERT1(dl.size() == dlsize || dl.size() == dlsize+1,
5265 "Wrong number of data for pointwize constraints brick");
5268 const model_complex_plain_vector *COEFF = 0;
5270 COEFF = &(md.complex_variable(dl[0]));
5272 GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
5273 "Data for coefficient should be a scalar");
5276 const model_complex_plain_vector &PT = md.complex_variable(dl[ind_pt]);
5277 size_type nb_co = gmm::vect_size(PT) / N;
5279 dim_type ind_unitv = dim_type((Q > 1) ? ind_pt+1 : 0);
5280 const model_complex_plain_vector &unitv
5281 = md.complex_variable(dl[ind_unitv]);
5282 GMM_ASSERT1((!ind_unitv || gmm::vect_size(unitv) == nb_co * Q),
5283 "Wrong size for vector of unit vectors");
5285 dim_type ind_rhs = dim_type((Q > 1) ? ind_pt+2 : ind_pt+1);
5286 if (dl.size() <
size_type(ind_rhs + 1)) ind_rhs = 0;
5287 const model_complex_plain_vector &rhs
5288 = md.complex_variable(dl[ind_rhs]);
5289 GMM_ASSERT1((!ind_rhs || gmm::vect_size(rhs) == nb_co),
5290 "Wrong size for vector of rhs");
5292 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
5293 || (penalized && (md.is_var_newer_than_brick(dl[ind_pt], ib)
5294 || md.is_var_newer_than_brick(dl[ind_unitv], ib)
5295 || md.is_var_newer_than_brick(dl[ind_rhs], ib)));
5297 if (recompute_matrix) {
5298 gmm::row_matrix<model_complex_sparse_vector>
5299 BB(nb_co*Q,mf_u.nb_dof());
5301 dal::bit_vector dof_untouched;
5302 getfem::mesh_trans_inv mti(mf_u.linked_mesh());
5306 (gmm::sub_vector(PT, gmm::sub_interval(i*N, N))), pt);
5309 gmm::row_matrix<model_complex_sparse_vector> &BBB = ((Q > 1) ? BB :cB);
5310 model_complex_plain_vector vv;
5311 interpolation(mf_u, mti, vv, vv, BBB, 1, 1, &dof_untouched);
5312 GMM_ASSERT1(dof_untouched.card() == 0,
5313 "Pointwize constraints : some of the points are outside "
5314 "the mesh: " << dof_untouched);
5319 gmm::add(gmm::scaled(gmm::mat_row(BB, i*Q+q), unitv[i*Q+q]),
5320 gmm::mat_row(cB, i));
5324 gmm::mult(gmm::transposed(cB), cB, matl[0]);
5325 gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
5333 gmm::mult(gmm::transposed(cB), rhs, vecl[0]);
5334 gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
5342 virtual std::string declare_volume_assembly_string
5343 (
const model &,
size_type,
const model::varnamelist &,
5344 const model::varnamelist &)
const {
5345 return std::string();
5348 pointwise_constraints_brick(
bool penalized) {
5349 set_flags(penalized ?
"Pointwise cosntraints with penalization brick"
5350 :
"Pointwise cosntraints with multipliers brick",
5360 (
model &md,
const std::string &varname,
5361 scalar_type penalisation_coeff,
const std::string &dataname_pt,
5362 const std::string &dataname_unitv,
const std::string &dataname_val) {
5363 std::string coeffname = md.
new_name(
"penalization_on_" + varname);
5369 pbrick pbr = std::make_shared<pointwise_constraints_brick>(
true);
5371 tl.push_back(model::term_description(varname, varname,
true));
5372 model::varnamelist vl(1, varname);
5373 model::varnamelist dl(1, coeffname);
5374 dl.push_back(dataname_pt);
5376 if (mf_u.
get_qdim() > 1) dl.push_back(dataname_unitv);
5377 if (dataname_val.size() > 0) dl.push_back(dataname_val);
5382 (
model &md,
const std::string &varname,
5383 const std::string &multname,
const std::string &dataname_pt,
5384 const std::string &dataname_unitv,
const std::string &dataname_val) {
5385 pbrick pbr = std::make_shared<pointwise_constraints_brick>(
false);
5387 tl.push_back(model::term_description(multname, varname,
true));
5388 model::varnamelist vl(1, varname);
5389 vl.push_back(multname);
5390 model::varnamelist dl(1, dataname_pt);
5392 if (mf_u.
get_qdim() > 1) dl.push_back(dataname_unitv);
5393 if (dataname_val.size() > 0) dl.push_back(dataname_val);
5398 (
model &md,
const std::string &varname,
const std::string &dataname_pt,
5399 const std::string &dataname_unitv,
const std::string &dataname_val) {
5400 std::string multname = md.
new_name(
"mult_on_" + varname);
5408 (md, varname, multname, dataname_pt, dataname_unitv, dataname_val);
5418 struct Helmholtz_brick :
public virtual_brick {
5420 virtual void asm_real_tangent_terms(
const model &md,
size_type,
5421 const model::varnamelist &vl,
5422 const model::varnamelist &dl,
5423 const model::mimlist &mims,
5424 model::real_matlist &matl,
5425 model::real_veclist &,
5426 model::real_veclist &,
5428 build_version)
const {
5429 GMM_ASSERT1(matl.size() == 1,
5430 "Helmholtz brick has one and only one term");
5431 GMM_ASSERT1(mims.size() == 1,
5432 "Helmholtz brick need one and only one mesh_im");
5433 GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
5434 "Wrong number of variables for Helmholtz brick");
5436 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5437 const mesh &m = mf_u.linked_mesh();
5439 GMM_ASSERT1(Q == 1,
"Helmholtz brick is only for scalar field, sorry.");
5440 const mesh_im &mim = *mims[0];
5441 const mesh_fem *mf_a = 0;
5442 mesh_region rg(region);
5443 m.intersect_with_mpi_region(rg);
5444 const model_real_plain_vector *A = &(md.real_variable(dl[0]));
5445 mf_a = md.pmesh_fem_of_variable(dl[0]);
5446 s = gmm::vect_size(*A);
5447 if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
5450 GMM_TRACE2(
"Stiffness matrix assembly for Helmholtz problem");
5451 gmm::clear(matl[0]);
5452 model_real_plain_vector A2(gmm::vect_size(*A));
5453 for (
size_type i=0; i < gmm::vect_size(*A); ++i)
5454 A2[i] = gmm::sqr((*A)[i]);
5460 GMM_ASSERT1(
false,
"Bad format Helmholtz brick coefficient");
5463 virtual void asm_complex_tangent_terms(
const model &md,
size_type,
5464 const model::varnamelist &vl,
5465 const model::varnamelist &dl,
5466 const model::mimlist &mims,
5467 model::complex_matlist &matl,
5468 model::complex_veclist &,
5469 model::complex_veclist &,
5471 build_version)
const {
5472 GMM_ASSERT1(matl.size() == 1,
5473 "Helmholtz brick has one and only one term");
5474 GMM_ASSERT1(mims.size() == 1,
5475 "Helmholtz brick need one and only one mesh_im");
5476 GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
5477 "Wrong number of variables for Helmholtz brick");
5479 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5480 const mesh &m = mf_u.linked_mesh();
5482 GMM_ASSERT1(Q == 1,
"Helmholtz brick is only for scalar field, sorry.");
5483 const mesh_im &mim = *mims[0];
5484 const mesh_fem *mf_a = 0;
5485 mesh_region rg(region);
5486 m.intersect_with_mpi_region(rg);
5487 const model_complex_plain_vector *A = &(md.complex_variable(dl[0]));
5488 mf_a = md.pmesh_fem_of_variable(dl[0]);
5489 s = gmm::vect_size(*A);
5490 if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
5493 GMM_TRACE2(
"Stiffness matrix assembly for Helmholtz problem");
5495 model_complex_plain_vector A2(gmm::vect_size(*A));
5496 for (
size_type i=0; i < gmm::vect_size(*A); ++i)
5497 A2[i] = gmm::sqr((*A)[i]);
5503 GMM_ASSERT1(
false,
"Bad format Helmholtz brick coefficient");
5507 set_flags(
"Helmholtz",
true ,
5515 const std::string &varname,
5516 const std::string &dataexpr,
5519 pbrick pbr = std::make_shared<Helmholtz_brick>();
5521 tl.push_back(model::term_description(varname, varname,
true));
5522 return md.
add_brick(pbr, model::varnamelist(1, varname),
5523 model::varnamelist(1, dataexpr), tl,
5524 model::mimlist(1, &mim), region);
5526 std::string test_varname
5527 =
"Test_" + sup_previous_and_dot_to_varname(varname);
5528 std::string expr =
"Grad_"+varname+
".Grad_"+test_varname
5529 +
" + sqr("+dataexpr+
")*"+varname+
"*"+test_varname;
5535 "Helmholtz (nonlinear)");
5548 struct Fourier_Robin_brick :
public virtual_brick {
5550 virtual void asm_real_tangent_terms(
const model &md,
size_type,
5551 const model::varnamelist &vl,
5552 const model::varnamelist &dl,
5553 const model::mimlist &mims,
5554 model::real_matlist &matl,
5555 model::real_veclist &,
5556 model::real_veclist &,
5558 build_version)
const {
5559 GMM_ASSERT1(matl.size() == 1,
5560 "Fourier-Robin brick has one and only one term");
5561 GMM_ASSERT1(mims.size() == 1,
5562 "Fourier-Robin brick need one and only one mesh_im");
5563 GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
5564 "Wrong number of variables for Fourier-Robin brick");
5566 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5567 const mesh &m = mf_u.linked_mesh();
5569 const mesh_im &mim = *mims[0];
5570 const mesh_fem *mf_a = 0;
5571 mesh_region rg(region);
5572 m.intersect_with_mpi_region(rg);
5573 const model_real_plain_vector *A = &(md.real_variable(dl[0]));
5574 mf_a = md.pmesh_fem_of_variable(dl[0]);
5575 s = gmm::vect_size(*A);
5576 if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
5577 GMM_ASSERT1(s == Q*Q,
5578 "Bad format Fourier-Robin brick coefficient");
5580 GMM_TRACE2(
"Fourier-Robin term assembly");
5581 gmm::clear(matl[0]);
5585 asm_homogeneous_qu_term(matl[0], mim, mf_u, *A, rg);
5588 virtual void asm_complex_tangent_terms(
const model &md,
size_type,
5589 const model::varnamelist &vl,
5590 const model::varnamelist &dl,
5591 const model::mimlist &mims,
5592 model::complex_matlist &matl,
5593 model::complex_veclist &,
5594 model::complex_veclist &,
5596 build_version)
const {
5597 GMM_ASSERT1(matl.size() == 1,
5598 "Fourier-Robin brick has one and only one term");
5599 GMM_ASSERT1(mims.size() == 1,
5600 "Fourier-Robin brick need one and only one mesh_im");
5601 GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
5602 "Wrong number of variables for Fourier-Robin brick");
5604 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5605 const mesh &m = mf_u.linked_mesh();
5607 const mesh_im &mim = *mims[0];
5608 const mesh_fem *mf_a = 0;
5609 mesh_region rg(region);
5610 m.intersect_with_mpi_region(rg);
5611 const model_complex_plain_vector *A = &(md.complex_variable(dl[0]));
5612 mf_a = md.pmesh_fem_of_variable(dl[0]);
5613 s = gmm::vect_size(*A);
5614 if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
5615 GMM_ASSERT1(s == Q*Q,
5616 "Bad format Fourier-Robin brick coefficient");
5618 GMM_TRACE2(
"Fourier-Robin term assembly");
5623 asm_homogeneous_qu_term(matl[0], mim, mf_u, *A, rg);
5626 Fourier_Robin_brick() {
5627 set_flags(
"Fourier Robin condition",
true ,
5636 const std::string &varname,
5637 const std::string &dataexpr,
5640 pbrick pbr = std::make_shared<Fourier_Robin_brick>();
5642 tl.push_back(model::term_description(varname, varname,
true));
5643 return md.
add_brick(pbr, model::varnamelist(1, varname),
5644 model::varnamelist(1, dataexpr), tl,
5645 model::mimlist(1, &mim), region);
5647 std::string test_varname
5648 =
"Test_" + sup_previous_and_dot_to_varname(varname);
5649 std::string expr =
"(("+dataexpr+
")*"+varname+
")."+test_varname;
5651 "Fourier-Robin",
true);
5654 "Fourier-Robin (nonlinear)");
5665 struct have_private_data_brick :
public virtual_brick {
5667 model_real_sparse_matrix rB;
5668 model_complex_sparse_matrix cB;
5669 model_real_plain_vector rL;
5670 model_complex_plain_vector cL;
5674 struct constraint_brick :
public have_private_data_brick {
5676 virtual void real_pre_assembly_in_serial(
const model &md,
size_type,
5677 const model::varnamelist &vl,
5678 const model::varnamelist &dl,
5679 const model::mimlist &mims,
5680 model::real_matlist &matl,
5681 model::real_veclist &vecl,
5682 model::real_veclist &,
5684 if (MPI_IS_MASTER()) {
5686 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5687 "Constraint brick has one and only one term");
5688 GMM_ASSERT1(mims.size() == 0,
5689 "Constraint brick need no mesh_im");
5690 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() <= 1,
5691 "Wrong number of variables for constraint brick");
5693 bool penalized = (vl.size() == 1);
5694 const model_real_plain_vector *COEFF = 0;
5696 bool has_data = (nameL.compare(
"") != 0);
5698 GMM_ASSERT1(nameL.compare(dl.back()) == 0 &&
5699 md.variable_exists(nameL) && md.is_data(nameL),
5701 const model_real_plain_vector &
5702 rrL = has_data ? md.real_variable(nameL) : rL;
5705 COEFF = &(md.real_variable(dl[0]));
5706 GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
5707 "Data for coefficient should be a scalar");
5710 gmm::scaled(rrL, gmm::abs((*COEFF)[0])), vecl[0]);
5712 gmm::scaled(rB, gmm::abs((*COEFF)[0])), matl[0]);
5720 virtual void complex_pre_assembly_in_serial(
const model &md,
size_type,
5721 const model::varnamelist &vl,
5722 const model::varnamelist &dl,
5723 const model::mimlist &mims,
5724 model::complex_matlist &matl,
5725 model::complex_veclist &vecl,
5726 model::complex_veclist &,
5728 build_version)
const {
5729 if (MPI_IS_MASTER()) {
5731 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5732 "Constraint brick has one and only one term");
5733 GMM_ASSERT1(mims.size() == 0,
5734 "Constraint brick need no mesh_im");
5735 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() <= 1,
5736 "Wrong number of variables for constraint brick");
5738 bool penalized = (vl.size() == 1);
5739 const model_complex_plain_vector *COEFF = 0;
5741 bool has_data = (nameL.compare(
"") != 0);
5743 GMM_ASSERT1(nameL.compare(dl.back()) == 0 &&
5744 md.variable_exists(nameL) && md.is_data(nameL),
5746 const model_complex_plain_vector &
5747 ccL = has_data ? md.complex_variable(nameL) : cL;
5750 COEFF = &(md.complex_variable(dl[0]));
5751 GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
5752 "Data for coefficient should be a scalar");
5755 gmm::scaled(ccL, gmm::abs((*COEFF)[0])), vecl[0]);
5757 gmm::scaled(cB, gmm::abs((*COEFF)[0])), matl[0]);
5765 virtual std::string declare_volume_assembly_string
5766 (
const model &,
size_type,
const model::varnamelist &,
5767 const model::varnamelist &)
const {
5768 return std::string();
5771 constraint_brick(
bool penalized) {
5772 set_flags(penalized ?
"Constraint with penalization brick"
5773 :
"Constraint with multipliers brick",
5782 model_real_sparse_matrix &set_private_data_brick_real_matrix
5784 pbrick pbr = md.brick_pointer(indbrick);
5785 md.touch_brick(indbrick);
5786 have_private_data_brick *p =
dynamic_cast<have_private_data_brick *
>
5787 (
const_cast<virtual_brick *
>(pbr.get()));
5788 GMM_ASSERT1(p,
"Wrong type of brick");
5792 model_real_plain_vector &set_private_data_brick_real_rhs
5794 pbrick pbr = md.brick_pointer(indbrick);
5795 md.touch_brick(indbrick);
5796 have_private_data_brick *p =
dynamic_cast<have_private_data_brick *
>
5797 (
const_cast<virtual_brick *
>(pbr.get()));
5798 GMM_ASSERT1(p,
"Wrong type of brick");
5799 if (p->nameL.compare(
"") != 0) GMM_WARNING1(
"Rhs already set by data name");
5803 model_complex_sparse_matrix &set_private_data_brick_complex_matrix
5805 pbrick pbr = md.brick_pointer(indbrick);
5806 md.touch_brick(indbrick);
5807 have_private_data_brick *p =
dynamic_cast<have_private_data_brick *
>
5808 (
const_cast<virtual_brick *
>(pbr.get()));
5809 GMM_ASSERT1(p,
"Wrong type of brick");
5813 model_complex_plain_vector &set_private_data_brick_complex_rhs
5815 pbrick pbr = md.brick_pointer(indbrick);
5816 md.touch_brick(indbrick);
5817 have_private_data_brick *p =
dynamic_cast<have_private_data_brick *
>
5818 (
const_cast<virtual_brick *
>(pbr.get()));
5819 GMM_ASSERT1(p,
"Wrong type of brick");
5820 if (p->nameL.compare(
"") != 0) GMM_WARNING1(
"Rhs already set by data name");
5824 void set_private_data_rhs
5825 (model &md,
size_type indbrick,
const std::string &varname) {
5826 pbrick pbr = md.brick_pointer(indbrick);
5827 md.touch_brick(indbrick);
5828 have_private_data_brick *p =
dynamic_cast<have_private_data_brick *
>
5829 (
const_cast<virtual_brick *
>(pbr.get()));
5830 GMM_ASSERT1(p,
"Wrong type of brick");
5831 if (p->nameL.compare(varname) != 0) {
5832 model::varnamelist dl = md.datanamelist_of_brick(indbrick);
5833 if (p->nameL.compare(
"") == 0) dl.push_back(varname);
5834 else dl.back() = varname;
5835 md.change_data_of_brick(indbrick, dl);
5840 size_type add_constraint_with_penalization
5841 (model &md,
const std::string &varname, scalar_type penalisation_coeff) {
5842 std::string coeffname = md.new_name(
"penalization_on_" + varname);
5843 md.add_fixed_size_data(coeffname, 1);
5844 if (md.is_complex())
5845 md.set_complex_variable(coeffname)[0] = penalisation_coeff;
5847 md.set_real_variable(coeffname)[0] = penalisation_coeff;
5848 pbrick pbr = std::make_shared<constraint_brick>(
true);
5850 tl.push_back(model::term_description(varname, varname,
true));
5851 model::varnamelist vl(1, varname);
5852 model::varnamelist dl(1, coeffname);
5853 return md.add_brick(pbr, vl, dl, tl, model::mimlist(),
size_type(-1));
5856 size_type add_constraint_with_multipliers
5857 (model &md,
const std::string &varname,
const std::string &multname) {
5858 pbrick pbr = std::make_shared<constraint_brick>(
false);
5860 tl.push_back(model::term_description(multname, varname,
true));
5861 model::varnamelist vl(1, varname);
5862 vl.push_back(multname);
5863 model::varnamelist dl;
5864 return md.add_brick(pbr, vl, dl, tl, model::mimlist(),
size_type(-1));
5874 struct explicit_matrix_brick :
public have_private_data_brick {
5876 virtual void real_pre_assembly_in_serial(
const model &,
size_type,
5877 const model::varnamelist &vl,
5878 const model::varnamelist &dl,
5879 const model::mimlist &mims,
5880 model::real_matlist &matl,
5881 model::real_veclist &vecl,
5882 model::real_veclist &,
5884 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5885 "Explicit matrix has one and only one term");
5886 GMM_ASSERT1(mims.size() == 0,
"Explicit matrix need no mesh_im");
5887 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() == 0,
5888 "Wrong number of variables for explicit matrix brick");
5889 GMM_ASSERT1(gmm::mat_ncols(rB) == gmm::mat_ncols(matl[0]) &&
5890 gmm::mat_nrows(rB) == gmm::mat_nrows(matl[0]),
5891 "Explicit matrix brick dimension mismatch ("<<
5892 gmm::mat_ncols(rB)<<
"x"<<gmm::mat_nrows(rB)<<
") != ("<<
5893 gmm::mat_ncols(matl[0])<<
"x"<<gmm::mat_nrows(matl[0])<<
")");
5897 virtual void complex_pre_assembly_in_serial(
const model &,
size_type,
5898 const model::varnamelist &vl,
5899 const model::varnamelist &dl,
5900 const model::mimlist &mims,
5901 model::complex_matlist &matl,
5902 model::complex_veclist &vecl,
5903 model::complex_veclist &,
5905 build_version)
const {
5906 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5907 "Explicit matrix has one and only one term");
5908 GMM_ASSERT1(mims.size() == 0,
"Explicit matrix need no mesh_im");
5909 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() == 0,
5910 "Wrong number of variables for explicit matrix brick");
5914 virtual std::string declare_volume_assembly_string
5915 (
const model &,
size_type,
const model::varnamelist &,
5916 const model::varnamelist &)
const {
5917 return std::string();
5920 explicit_matrix_brick(
bool symmetric_,
bool coercive_) {
5921 set_flags(
"Explicit matrix brick",
5923 symmetric_ , coercive_ ,
5930 (model &md,
const std::string &varname1,
const std::string &varname2,
5931 bool issymmetric,
bool iscoercive) {
5932 pbrick pbr = std::make_shared<explicit_matrix_brick>(issymmetric,
5935 tl.push_back(model::term_description(varname1, varname2, issymmetric));
5936 model::varnamelist vl(1, varname1);
5937 vl.push_back(varname2);
5938 model::varnamelist dl;
5939 return md.add_brick(pbr, vl, dl, tl, model::mimlist(),
size_type(-1));
5948 struct explicit_rhs_brick :
public have_private_data_brick {
5950 virtual void real_pre_assembly_in_serial(
const model &,
size_type,
5951 const model::varnamelist &vl,
5952 const model::varnamelist &dl,
5953 const model::mimlist &mims,
5954 model::real_matlist &matl,
5955 model::real_veclist &vecl,
5956 model::real_veclist &,
5958 if (MPI_IS_MASTER()) {
5959 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5960 "Explicit rhs has one and only one term");
5961 GMM_ASSERT1(mims.size() == 0,
"Explicit rhs need no mesh_im");
5962 GMM_ASSERT1(vl.size() == 1 && dl.size() == 0,
5963 "Wrong number of variables for explicit rhs brick");
5968 virtual void complex_pre_assembly_in_serial(
const model &,
size_type,
5969 const model::varnamelist &vl,
5970 const model::varnamelist &dl,
5971 const model::mimlist &mims,
5972 model::complex_matlist &matl,
5973 model::complex_veclist &vecl,
5974 model::complex_veclist &,
5976 build_version)
const {
5977 if (MPI_IS_MASTER()) {
5978 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5979 "Explicit rhs has one and only one term");
5980 GMM_ASSERT1(mims.size() == 0,
"Explicit rhs need no mesh_im");
5981 GMM_ASSERT1(vl.size() == 1 && dl.size() == 0,
5982 "Wrong number of variables for explicit rhs brick");
5988 virtual std::string declare_volume_assembly_string
5989 (
const model &,
size_type,
const model::varnamelist &,
5990 const model::varnamelist &)
const {
5991 return std::string();
5994 explicit_rhs_brick() {
5995 set_flags(
"Explicit rhs brick",
6005 (model &md,
const std::string &varname) {
6006 pbrick pbr = std::make_shared<explicit_rhs_brick>();
6008 tl.push_back(model::term_description(varname));
6009 model::varnamelist vl(1, varname);
6010 model::varnamelist dl;
6011 return md.add_brick(pbr, vl, dl, tl, model::mimlist(),
size_type(-1));
6021 struct iso_lin_elasticity_new_brick :
public virtual_brick {
6023 std::string expr, dataname3;
6025 void asm_real_tangent_terms(
const model &md,
size_type ib,
6026 const model::varnamelist &vl,
6027 const model::varnamelist &dl,
6028 const model::mimlist &mims,
6029 model::real_matlist &matl,
6030 model::real_veclist &vecl,
6031 model::real_veclist &,
6033 build_version version)
const override {
6034 GMM_ASSERT1(vl.size() == 1,
"Linearized isotropic elasticity brick "
6035 "has one and only one variable");
6036 GMM_ASSERT1(matl.size() == 1,
"Linearized isotropic elasticity brick "
6037 "has one and only one term");
6038 GMM_ASSERT1(mims.size() == 1,
"Linearized isotropic elasticity brick "
6039 "needs one and only one mesh_im");
6041 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0);
6042 for (
size_type i = 0; i < dl.size(); ++i) {
6043 recompute_matrix = recompute_matrix ||
6044 md.is_var_newer_than_brick(dl[i], ib);
6047 if (recompute_matrix) {
6049 ga_workspace workspace(md, ga_workspace::inherit::ALL);
6050 workspace.add_expression(expr, *(mims[0]), region);
6051 GMM_TRACE2(name <<
": generic matrix assembly");
6052 workspace.assembly(2);
6053 scalar_type
alpha = scalar_type(1)
6054 / (workspace.factor_of_variable(vl[0]));
6055 const auto &R=workspace.assembled_matrix();
6056 gmm::sub_interval I = workspace.interval_of_variable(vl[0]);
6057 gmm::copy(gmm::scaled(gmm::sub_matrix(R, I, I), alpha),
6061 if (dataname3.size()) {
6067 gmm::scaled(md.real_variable(dataname3), scalar_type(-1)),
6073 void real_post_assembly_in_serial(
const model &md,
size_type ib,
6074 const model::varnamelist &,
6075 const model::varnamelist &,
6076 const model::mimlist &,
6077 model::real_matlist &,
6078 model::real_veclist &vecl,
6079 model::real_veclist &,
6081 build_version)
const override {
6082 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
6086 virtual std::string declare_volume_assembly_string
6087 (
const model &,
size_type,
const model::varnamelist &,
6088 const model::varnamelist &)
const {
6092 iso_lin_elasticity_new_brick(
const std::string &expr_,
6093 const std::string &dataname3_) {
6094 expr = expr_; dataname3 = dataname3_;
6095 set_flags(
"Linearized isotropic elasticity",
true ,
6104 (
model &md,
const mesh_im &mim,
const std::string &varname,
6105 const std::string &dataexpr1,
const std::string &dataexpr2,
6106 size_type region,
const std::string &dataname3) {
6107 std::string test_varname
6108 =
"Test_" + sup_previous_and_dot_to_varname(varname);
6110 std::string expr1 =
"((("+dataexpr1+
")*(Div_"+varname+
"-Div_"+dataname3
6111 +
"))*Id(meshdim)+(2*("+dataexpr2+
"))*(Sym(Grad_"+varname
6112 +
")-Sym(Grad_"+dataname3+
"))):Grad_" +test_varname;
6113 std::string expr2 =
"(Div_"+varname+
"*(("+dataexpr1+
")*Id(meshdim))"
6114 +
"+(2*("+dataexpr2+
"))*Sym(Grad_"+varname+
")):Grad_"+test_varname;
6117 model::varnamelist vl, dl;
6119 ga_workspace workspace(md, ga_workspace::inherit::ALL);
6120 workspace.add_expression(expr2, mim, region);
6121 model::varnamelist vl_test1, vl_test2;
6122 is_lin = workspace.used_variables(vl, vl_test1, vl_test2, dl, 2);
6125 pbrick pbr = std::make_shared<iso_lin_elasticity_new_brick>
6128 tl.push_back(model::term_description(varname,
6129 sup_previous_and_dot_to_varname(varname),
true));
6130 if (dataname3.size()) dl.push_back(dataname3);
6131 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
6134 (md, mim, dataname3.size() ? expr1 : expr2, region,
false,
false,
6135 "Linearized isotropic elasticity (with nonlinear dependance)");
6140 (
model &md,
const mesh_im &mim,
const std::string &varname,
6141 const std::string &data_E,
const std::string &data_nu,
6143 std::string test_varname
6144 =
"Test_" + sup_previous_and_dot_to_varname(varname);
6146 std::string mu =
"(("+data_E+
")/(2*(1+("+data_nu+
"))))";
6147 std::string lambda =
"(("+data_E+
")*("+data_nu+
")/((1+("+data_nu+
"))*(1-2*("
6149 std::string expr = lambda+
"*Div_"+varname+
"*Div_"+test_varname
6150 +
"+"+mu+
"*(Grad_"+varname+
"+Grad_"+varname+
"'):Grad_"+test_varname;
6154 ga_workspace workspace(md, ga_workspace::inherit::ALL);
6155 workspace.add_expression(expr, mim, region);
6156 is_lin = workspace.is_linear(2);
6160 "Linearized isotropic elasticity");
6163 (md, mim, expr, region,
false,
false,
6164 "Linearized isotropic elasticity (with nonlinear dependance)");
6169 (
model &md,
const mesh_im &mim,
const std::string &varname,
6170 const std::string &data_E,
const std::string &data_nu,
6172 std::string test_varname
6173 =
"Test_" + sup_previous_and_dot_to_varname(varname);
6176 GMM_ASSERT1(mfu,
"The variable should be a fem variable");
6179 std::string mu =
"(("+data_E+
")/(2*(1+("+data_nu+
"))))";
6180 std::string lambda =
"(("+data_E+
")*("+data_nu+
")/((1+("+data_nu+
"))*(1-2*("
6183 lambda =
"(("+data_E+
")*("+data_nu+
")/((1-sqr("+data_nu+
"))))";
6184 std::string expr = lambda+
"*Div_"+varname+
"*Div_"+test_varname
6185 +
"+"+mu+
"*(Grad_"+varname+
"+Grad_"+varname+
"'):Grad_"+test_varname;
6189 ga_workspace workspace(md, ga_workspace::inherit::ALL);
6190 workspace.add_expression(expr, mim, region);
6191 is_lin = workspace.is_linear(2);
6195 "Linearized isotropic elasticity");
6198 (md, mim, expr, region,
false,
false,
6199 "Linearized isotropic elasticity (with nonlinear dependance)");
6204 void compute_isotropic_linearized_Von_Mises_or_Tresca
6205 (model &md,
const std::string &varname,
const std::string &data_lambda,
6206 const std::string &data_mu,
const mesh_fem &mf_vm,
6207 model_real_plain_vector &VM,
bool tresca) {
6210 const mesh_fem &mf_u = md.mesh_fem_of_variable(varname);
6211 const mesh_fem *mf_lambda = md.pmesh_fem_of_variable(data_lambda);
6212 const model_real_plain_vector *lambda=&(md.real_variable(data_lambda));
6213 const mesh_fem *mf_mu = md.pmesh_fem_of_variable(data_mu);
6214 const model_real_plain_vector *mu = &(md.real_variable(data_mu));
6217 if (mf_lambda) sl = sl * mf_lambda->get_qdim() / mf_lambda->nb_dof();
6219 if (mf_mu) sm = sm * mf_mu->get_qdim() / mf_mu->nb_dof();
6221 GMM_ASSERT1(sl == 1 && sm == 1,
"Bad format for Lame coefficients");
6222 GMM_ASSERT1(mf_lambda == mf_mu,
6223 "The two Lame coefficients should be described on the same "
6224 "finite element method.");
6228 md.real_variable(varname), VM,
6229 *mf_lambda, *lambda,
6234 model_real_plain_vector LAMBDA(mf_lambda->nb_dof(), (*lambda)[0]);
6235 model_real_plain_vector MU(mf_lambda->nb_dof(), (*mu)[0]);
6237 md.real_variable(varname), VM,
6246 std::string sigma_d=
"("+data_mu+
")*(Grad_"+varname+
"+Grad_"+varname+
"')";
6247 std::string expr =
"sqrt(3/2)*Norm(Deviator("+sigma_d+
"))";
6248 ga_interpolation_Lagrange_fem(md, expr, mf_vm, VM);
6254 (
model &md,
const std::string &varname,
const std::string &data_E,
6255 const std::string &data_nu,
const mesh_fem &mf_vm,
6256 model_real_plain_vector &VM) {
6260 std::string mu =
"(("+data_E+
")/(2*(1+("+data_nu+
"))))";
6263 std::string sigma_d = mu+
"*(Grad_"+varname+
"+Grad_"+varname+
"')";
6264 std::string expr =
"sqrt(3/2)*Norm(Deviator("+sigma_d+
"))";
6265 ga_interpolation_Lagrange_fem(md, expr, mf_vm, VM);
6269 (
model &md,
const std::string &varname,
const std::string &data_E,
6270 const std::string &data_nu,
const mesh_fem &mf_vm,
6271 model_real_plain_vector &VM) {
6280 std::string mu =
"(("+data_E+
")/(2*(1+("+data_nu+
"))))";
6283 std::string sigma_d = mu+
"*(Grad_"+varname+
"+Grad_"+varname+
"')";
6284 std::string expr =
"sqrt(3/2)*Norm(Deviator("+sigma_d+
"))";
6285 ga_interpolation_Lagrange_fem(md, expr, mf_vm, VM);
6295 struct linear_incompressibility_brick :
public virtual_brick {
6297 virtual void asm_real_tangent_terms(
const model &md,
size_type ,
6298 const model::varnamelist &vl,
6299 const model::varnamelist &dl,
6300 const model::mimlist &mims,
6301 model::real_matlist &matl,
6302 model::real_veclist &,
6303 model::real_veclist &,
6305 build_version)
const {
6307 GMM_ASSERT1((matl.size() == 1 && dl.size() == 0)
6308 || (matl.size() == 2 && dl.size() == 1),
6309 "Wrong term and/or data number for Linear incompressibility "
6311 GMM_ASSERT1(mims.size() == 1,
"Linear incompressibility brick need one "
6312 "and only one mesh_im");
6313 GMM_ASSERT1(vl.size() == 2,
"Wrong number of variables for linear "
6314 "incompressibility brick");
6316 bool penalized = (dl.size() == 1);
6317 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6318 const mesh_fem &mf_p = md.mesh_fem_of_variable(vl[1]);
6319 const mesh_im &mim = *mims[0];
6320 const model_real_plain_vector *COEFF = 0;
6321 const mesh_fem *mf_data = 0;
6324 COEFF = &(md.real_variable(dl[0]));
6325 mf_data = md.pmesh_fem_of_variable(dl[0]);
6327 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
6328 GMM_ASSERT1(s == 1,
"Bad format for the penalization parameter");
6331 mesh_region rg(region);
6332 mim.linked_mesh().intersect_with_mpi_region(rg);
6334 GMM_TRACE2(
"Stokes term assembly");
6342 gmm::scale(matl[1], scalar_type(-1));
6346 gmm::scale(matl[1], -(*COEFF)[0]);
6353 virtual void real_post_assembly_in_serial(
const model &,
size_type,
6354 const model::varnamelist &,
6355 const model::varnamelist &,
6356 const model::mimlist &,
6357 model::real_matlist &,
6358 model::real_veclist &,
6359 model::real_veclist &,
6361 build_version)
const
6365 linear_incompressibility_brick() {
6366 set_flags(
"Linear incompressibility brick",
6375 (
model &md,
const mesh_im &mim,
const std::string &varname,
6376 const std::string &multname,
size_type region,
6377 const std::string &dataexpr) {
6379 pbrick pbr = std::make_shared<linear_incompressibility_brick>();
6381 tl.push_back(model::term_description(multname, varname,
true));
6382 model::varnamelist vl(1, varname);
6383 vl.push_back(multname);
6384 model::varnamelist dl;
6385 if (dataname.size()) {
6386 dl.push_back(dataexpr);
6387 tl.push_back(model::term_description(multname, multname,
true));
6389 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
6391 std::string test_varname
6392 =
"Test_" + sup_previous_and_dot_to_varname(varname);
6393 std::string test_multname
6394 =
"Test_" + sup_previous_and_dot_to_varname(multname);
6396 if (dataexpr.size())
6397 expr =
"-"+multname+
"*Div_"+test_varname +
"-"+test_multname
6398 +
"*Div_"+varname+
"+(("+dataexpr+
")*"+multname+
")*"+test_multname;
6400 expr =
"-"+multname+
"*Div_"+test_varname +
"-"+test_multname
6403 "Linear incompressibility",
true);
6406 (md, mim, expr, region,
false,
false,
6407 "Linear incompressibility (with nonlinear dependance)");
6420 struct mass_brick :
public virtual_brick {
6422 virtual void asm_real_tangent_terms(
const model &md,
size_type,
6423 const model::varnamelist &vl,
6424 const model::varnamelist &dl,
6425 const model::mimlist &mims,
6426 model::real_matlist &matl,
6427 model::real_veclist &,
6428 model::real_veclist &,
6430 build_version)
const {
6431 GMM_ASSERT1(matl.size() == 1,
6432 "Mass brick has one and only one term");
6433 GMM_ASSERT1(mims.size() == 1,
6434 "Mass brick need one and only one mesh_im");
6435 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
6436 "Wrong number of variables for mass brick");
6438 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6439 const mesh &m = mf_u.linked_mesh();
6440 const mesh_im &mim = *mims[0];
6441 mesh_region rg(region);
6442 m.intersect_with_mpi_region(rg);
6444 const mesh_fem *mf_rho = 0;
6445 const model_real_plain_vector *rho = 0;
6448 mf_rho = md.pmesh_fem_of_variable(dl[0]);
6449 rho = &(md.real_variable(dl[0]));
6451 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6452 GMM_ASSERT1(sl == 1,
"Bad format of mass brick coefficient");
6455 GMM_TRACE2(
"Mass matrix assembly");
6457 if (dl.size() && mf_rho) {
6461 if (dl.size()) gmm::scale(matl[0], (*rho)[0]);
6465 virtual void asm_complex_tangent_terms(
const model &md,
size_type,
6466 const model::varnamelist &vl,
6467 const model::varnamelist &dl,
6468 const model::mimlist &mims,
6469 model::complex_matlist &matl,
6470 model::complex_veclist &,
6471 model::complex_veclist &,
6473 build_version)
const {
6474 GMM_ASSERT1(matl.size() == 1,
6475 "Mass brick has one and only one term");
6476 GMM_ASSERT1(mims.size() == 1,
6477 "Mass brick need one and only one mesh_im");
6478 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
6479 "Wrong number of variables for mass brick");
6481 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6482 const mesh &m = mf_u.linked_mesh();
6483 const mesh_im &mim = *mims[0];
6484 mesh_region rg(region);
6485 m.intersect_with_mpi_region(rg);
6487 const mesh_fem *mf_rho = 0;
6488 const model_complex_plain_vector *rho = 0;
6491 mf_rho = md.pmesh_fem_of_variable(dl[0]);
6492 rho = &(md.complex_variable(dl[0]));
6494 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6495 GMM_ASSERT1(sl == 1,
"Bad format of mass brick coefficient");
6498 GMM_TRACE2(
"Mass matrix assembly");
6500 if (dl.size() && mf_rho) {
6504 if (dl.size()) gmm::scale(matl[0], (*rho)[0]);
6508 virtual std::string declare_volume_assembly_string
6509 (
const model &,
size_type,
const model::varnamelist &,
6510 const model::varnamelist &)
const {
6511 return std::string();
6515 set_flags(
"Mass brick",
true ,
6524 (
model &md,
const mesh_im &mim,
const std::string &varname,
6525 const std::string &dataexpr_rho,
size_type region) {
6527 pbrick pbr = std::make_shared<mass_brick>();
6529 tl.push_back(model::term_description(varname, varname,
true));
6530 model::varnamelist dl;
6531 if (dataexpr_rho.size())
6532 dl.push_back(dataexpr_rho);
6533 return md.
add_brick(pbr, model::varnamelist(1, varname), dl, tl,
6534 model::mimlist(1, &mim), region);
6536 std::string test_varname
6537 =
"Test_" + sup_previous_and_dot_to_varname(varname);
6539 if (dataexpr_rho.size())
6540 expr =
"(("+dataexpr_rho+
")*"+varname+
")."+test_varname;
6542 expr = varname+
"."+test_varname;
6544 "Mass matrix",
true);
6547 "Mass matrix (nonlinear)");
6558 struct lumped_mass_for_first_order_brick :
public virtual_brick {
6560 virtual void asm_real_tangent_terms(
const model &md,
size_type,
6561 const model::varnamelist &vl,
6562 const model::varnamelist &dl,
6563 const model::mimlist &mims,
6564 model::real_matlist &matl,
6565 model::real_veclist &,
6566 model::real_veclist &,
6568 build_version)
const {
6569 GMM_ASSERT1(matl.size() == 1,
6570 "Lumped Mass brick has one and only one term");
6571 GMM_ASSERT1(mims.size() == 1,
6572 "Lumped Mass brick needs one and only one mesh_im");
6573 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
6574 "Wrong number of variables for lumped mass brick");
6576 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6577 const mesh &m = mf_u.linked_mesh();
6578 const mesh_im &mim = *mims[0];
6579 mesh_region rg(region);
6580 m.intersect_with_mpi_region(rg);
6582 const mesh_fem *mf_rho = 0;
6583 const model_real_plain_vector *rho = 0;
6586 mf_rho = md.pmesh_fem_of_variable(dl[0]);
6587 rho = &(md.real_variable(dl[0]));
6589 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6590 GMM_ASSERT1(sl == 1,
"Bad format of mass brick coefficient");
6593 GMM_TRACE2(
"Lumped mass matrix assembly (please check that integration is 1st order.)");
6595 if (dl.size() && mf_rho) {
6599 if (dl.size()) gmm::scale(matl[0], (*rho)[0]);
6604 lumped_mass_for_first_order_brick() {
6605 set_flags(
"Lumped mass brick",
true ,
6614 (
model & md,
const mesh_im &mim,
const std::string &varname,
6615 const std::string &dataexpr_rho,
size_type region) {
6616 pbrick pbr = std::make_shared<lumped_mass_for_first_order_brick>();
6618 tl.push_back(model::term_description(varname, varname,
true));
6619 model::varnamelist dl;
6620 if (dataexpr_rho.size())
6621 dl.push_back(dataexpr_rho);
6622 return md.
add_brick(pbr, model::varnamelist(1, varname), dl, tl,
6623 model::mimlist(1, &mim), region);
6639 struct basic_d_on_dt_brick :
public virtual_brick {
6641 virtual void asm_real_tangent_terms(
const model &md,
size_type ib,
6642 const model::varnamelist &vl,
6643 const model::varnamelist &dl,
6644 const model::mimlist &mims,
6645 model::real_matlist &matl,
6646 model::real_veclist &vecl,
6647 model::real_veclist &,
6649 build_version version)
const {
6650 GMM_ASSERT1(matl.size() == 1,
6651 "Basic d/dt brick has one and only one term");
6652 GMM_ASSERT1(mims.size() == 1,
6653 "Basic d/dt brick need one and only one mesh_im");
6654 GMM_ASSERT1(vl.size() == 1 && dl.size() >= 2 && dl.size() <= 3,
6655 "Wrong number of variables for basic d/dt brick");
6659 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
6660 || (md.is_var_newer_than_brick(dl[1], ib));
6662 recompute_matrix = recompute_matrix ||
6663 md.is_var_newer_than_brick(dl[2], ib);
6666 if (recompute_matrix) {
6667 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6668 const mesh &m = mf_u.linked_mesh();
6669 const mesh_im &mim = *mims[0];
6670 mesh_region rg(region);
6671 m.intersect_with_mpi_region(rg);
6673 const model_real_plain_vector &dt = md.real_variable(dl[1]);
6674 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for time step");
6676 const mesh_fem *mf_rho = 0;
6677 const model_real_plain_vector *rho = 0;
6679 if (dl.size() > 2) {
6680 mf_rho = md.pmesh_fem_of_variable(dl[2]);
6681 rho = &(md.real_variable(dl[2]));
6683 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6684 GMM_ASSERT1(sl == 1,
"Bad format for density");
6687 GMM_TRACE2(
"Mass matrix assembly for d_on_dt brick");
6688 if (dl.size() > 2 && mf_rho) {
6691 gmm::scale(matl[0], scalar_type(1) / dt[0]);
6695 if (dl.size() > 2) gmm::scale(matl[0], (*rho)[0] / dt[0]);
6696 else gmm::scale(matl[0], scalar_type(1) / dt[0]);
6699 gmm::mult(matl[0], md.real_variable(dl[0], 1), vecl[0]);
6702 virtual void asm_complex_tangent_terms(
const model &md,
size_type ib,
6703 const model::varnamelist &vl,
6704 const model::varnamelist &dl,
6705 const model::mimlist &mims,
6706 model::complex_matlist &matl,
6707 model::complex_veclist &vecl,
6708 model::complex_veclist &,
6710 build_version version)
const {
6711 GMM_ASSERT1(matl.size() == 1,
6712 "Basic d/dt brick has one and only one term");
6713 GMM_ASSERT1(mims.size() == 1,
6714 "Basic d/dt brick need one and only one mesh_im");
6715 GMM_ASSERT1(vl.size() == 1 && dl.size() >= 2 && dl.size() <= 3,
6716 "Wrong number of variables for basic d/dt brick");
6720 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
6721 || (md.is_var_newer_than_brick(dl[1], ib));
6723 recompute_matrix = recompute_matrix ||
6724 md.is_var_newer_than_brick(dl[2], ib);
6726 if (recompute_matrix) {
6727 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6728 const mesh &m = mf_u.linked_mesh();
6729 const mesh_im &mim = *mims[0];
6731 mesh_region rg(region);
6732 m.intersect_with_mpi_region(rg);
6734 const model_complex_plain_vector &dt = md.complex_variable(dl[1]);
6735 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for time step");
6737 const mesh_fem *mf_rho = 0;
6738 const model_complex_plain_vector *rho = 0;
6740 if (dl.size() > 2) {
6741 mf_rho = md.pmesh_fem_of_variable(dl[2]);
6742 rho = &(md.complex_variable(dl[2]));
6744 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6745 GMM_ASSERT1(sl == 1,
"Bad format for density");
6748 GMM_TRACE2(
"Mass matrix assembly for d_on_dt brick");
6749 if (dl.size() > 2 && mf_rho) {
6752 gmm::scale(matl[0], scalar_type(1) / dt[0]);
6756 if (dl.size() > 2) gmm::scale(matl[0], (*rho)[0] / dt[0]);
6757 else gmm::scale(matl[0], scalar_type(1) / dt[0]);
6760 gmm::mult(matl[0], md.complex_variable(dl[0], 1), vecl[0]);
6763 virtual std::string declare_volume_assembly_string
6764 (
const model &,
size_type,
const model::varnamelist &,
6765 const model::varnamelist &)
const {
6766 return std::string();
6769 basic_d_on_dt_brick() {
6770 set_flags(
"Basic d/dt brick",
true ,
6779 (
model &md,
const mesh_im &mim,
const std::string &varname,
6780 const std::string &dataname_dt,
const std::string &dataname_rho,
6782 pbrick pbr = std::make_shared<basic_d_on_dt_brick>();
6784 tl.push_back(model::term_description(varname, varname,
true));
6785 model::varnamelist dl(1, varname);
6786 dl.push_back(dataname_dt);
6787 if (dataname_rho.size())
6788 dl.push_back(dataname_rho);
6789 return md.
add_brick(pbr, model::varnamelist(1, varname), dl, tl,
6790 model::mimlist(1, &mim), region);
6801 struct basic_d2_on_dt2_brick :
public virtual_brick {
6803 mutable scalar_type old_alphadt2;
6805 virtual void asm_real_tangent_terms(
const model &md,
size_type ib,
6806 const model::varnamelist &vl,
6807 const model::varnamelist &dl,
6808 const model::mimlist &mims,
6809 model::real_matlist &matl,
6810 model::real_veclist &vecl,
6811 model::real_veclist &,
6813 build_version version)
const {
6814 GMM_ASSERT1(matl.size() == 1,
6815 "Basic d2/dt2 brick has one and only one term");
6816 GMM_ASSERT1(mims.size() == 1,
6817 "Basic d2/dt2 brick need one and only one mesh_im");
6818 GMM_ASSERT1(vl.size() == 1 && dl.size() >= 4 && dl.size() <= 5,
6819 "Wrong number of variables for basic d2/dt2 brick");
6821 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0);
6823 recompute_matrix = recompute_matrix || md.is_var_newer_than_brick(dl[2], ib);
6824 if (dl.size() > 4) recompute_matrix || md.is_var_newer_than_brick(dl[4], ib);
6826 const model_real_plain_vector &dt = md.real_variable(dl[2]);
6827 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for time step");
6828 const model_real_plain_vector &alpha = md.real_variable(dl[3]);
6829 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for parameter alpha");
6830 scalar_type alphadt2 = gmm::sqr(dt[0]) * alpha[0];
6832 if (!recompute_matrix && alphadt2 != old_alphadt2)
6833 gmm::scale(matl[0], old_alphadt2/alphadt2);
6834 old_alphadt2 = alphadt2;
6836 if (recompute_matrix) {
6837 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6838 const mesh &m = mf_u.linked_mesh();
6839 const mesh_im &mim = *mims[0];
6840 mesh_region rg(region);
6841 m.intersect_with_mpi_region(rg);
6843 const mesh_fem *mf_rho = 0;
6844 const model_real_plain_vector *rho = 0;
6846 if (dl.size() > 4) {
6847 mf_rho = md.pmesh_fem_of_variable(dl[4]);
6848 rho = &(md.real_variable(dl[4]));
6850 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6851 GMM_ASSERT1(sl == 1,
"Bad format for density");
6854 GMM_TRACE2(
"Mass matrix assembly for d2_on_dt2 brick");
6855 if (dl.size() > 4 && mf_rho) {
6858 gmm::scale(matl[0], scalar_type(1) / alphadt2);
6862 if (dl.size() > 4) gmm::scale(matl[0], (*rho)[0] / alphadt2);
6863 else gmm::scale(matl[0], scalar_type(1) / alphadt2);
6866 gmm::mult(matl[0], md.real_variable(dl[0], 1), vecl[0]);
6867 gmm::mult_add(matl[0], gmm::scaled(md.real_variable(dl[1], 1), dt[0]),
6871 virtual void asm_complex_tangent_terms(
const model &md,
size_type ib,
6872 const model::varnamelist &vl,
6873 const model::varnamelist &dl,
6874 const model::mimlist &mims,
6875 model::complex_matlist &matl,
6876 model::complex_veclist &vecl,
6877 model::complex_veclist &,
6879 build_version version)
const {
6880 GMM_ASSERT1(matl.size() == 1,
6881 "Basic d2/dt2 brick has one and only one term");
6882 GMM_ASSERT1(mims.size() == 1,
6883 "Basic d2/dt2 brick need one and only one mesh_im");
6884 GMM_ASSERT1(vl.size() == 1 && dl.size() >= 4 && dl.size() <= 5,
6885 "Wrong number of variables for basic d2/dt2 brick");
6887 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0);
6889 recompute_matrix = recompute_matrix || md.is_var_newer_than_brick(dl[2], ib);
6890 if (dl.size() > 4) recompute_matrix || md.is_var_newer_than_brick(dl[4], ib);
6893 const model_complex_plain_vector &dt = md.complex_variable(dl[2]);
6894 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for time step");
6895 const model_complex_plain_vector &
alpha = md.complex_variable(dl[3]);
6896 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for parameter alpha");
6897 scalar_type alphadt2 = gmm::real(gmm::sqr(dt[0]) * alpha[0]);
6899 if (!recompute_matrix && alphadt2 != old_alphadt2)
6900 gmm::scale(matl[0], old_alphadt2/alphadt2);
6901 old_alphadt2 = alphadt2;
6903 if (recompute_matrix) {
6904 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6905 const mesh &m = mf_u.linked_mesh();
6906 const mesh_im &mim = *mims[0];
6907 mesh_region rg(region);
6908 m.intersect_with_mpi_region(rg);
6910 const mesh_fem *mf_rho = 0;
6911 const model_complex_plain_vector *rho = 0;
6913 if (dl.size() > 4) {
6914 mf_rho = md.pmesh_fem_of_variable(dl[4]);
6915 rho = &(md.complex_variable(dl[4]));
6917 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6918 GMM_ASSERT1(sl == 1,
"Bad format for density");
6921 GMM_TRACE2(
"Mass matrix assembly for d2_on_dt2 brick");
6922 if (dl.size() > 4 && mf_rho) {
6925 gmm::scale(matl[0], scalar_type(1) / alphadt2);
6929 if (dl.size() > 4) gmm::scale(matl[0], (*rho)[0] / alphadt2);
6930 else gmm::scale(matl[0], scalar_type(1) / alphadt2);
6933 gmm::mult(matl[0], md.complex_variable(dl[0], 1), vecl[0]);
6934 gmm::mult_add(matl[0], gmm::scaled(md.complex_variable(dl[1], 1), dt[0]),
6938 virtual std::string declare_volume_assembly_string
6939 (
const model &,
size_type,
const model::varnamelist &,
6940 const model::varnamelist &)
const {
6941 return std::string();
6944 basic_d2_on_dt2_brick() {
6945 set_flags(
"Basic d2/dt2 brick",
true ,
6954 (
model &md,
const mesh_im &mim,
const std::string &varnameU,
6955 const std::string &datanameV,
6956 const std::string &dataname_dt,
6957 const std::string &dataname_alpha,
6958 const std::string &dataname_rho,
6960 pbrick pbr = std::make_shared<basic_d2_on_dt2_brick>();
6962 tl.push_back(model::term_description(varnameU, varnameU,
true));
6963 model::varnamelist dl(1, varnameU);
6964 dl.push_back(datanameV);
6965 dl.push_back(dataname_dt);
6966 dl.push_back(dataname_alpha);
6967 if (dataname_rho.size())
6968 dl.push_back(dataname_rho);
6969 return md.
add_brick(pbr, model::varnamelist(1, varnameU), dl, tl,
6970 model::mimlist(1, &mim), region);
6989 void theta_method_dispatcher::set_dispatch_coeff(
const model &md,
size_type ib)
const {
6991 if (md.is_complex())
6992 theta = gmm::real(md.complex_variable(param_names[0])[0]);
6994 theta = md.real_variable(param_names[0])[0];
6996 md.matrix_coeff_of_brick(ib) = theta;
6998 md.rhs_coeffs_of_brick(ib)[0] = theta;
7000 md.rhs_coeffs_of_brick(ib)[1] = (scalar_type(1) - theta);
7004 void theta_method_dispatcher::next_real_iter
7005 (
const model &md,
size_type ib,
const model::varnamelist &vl,
7006 const model::varnamelist &dl, model::real_matlist &matl,
7007 std::vector<model::real_veclist> &vectl,
7008 std::vector<model::real_veclist> &vectl_sym,
bool first_iter)
const {
7009 next_iter(md, ib, vl, dl, matl, vectl, vectl_sym, first_iter);
7012 void theta_method_dispatcher::next_complex_iter
7013 (
const model &md,
size_type ib,
const model::varnamelist &vl,
7014 const model::varnamelist &dl,
7015 model::complex_matlist &matl,
7016 std::vector<model::complex_veclist> &vectl,
7017 std::vector<model::complex_veclist> &vectl_sym,
7018 bool first_iter)
const {
7019 next_iter(md, ib, vl, dl, matl, vectl, vectl_sym, first_iter);
7022 void theta_method_dispatcher::asm_real_tangent_terms
7023 (
const model &md,
size_type ib, model::real_matlist &,
7024 std::vector<model::real_veclist> &,
7025 std::vector<model::real_veclist> &,
7026 build_version version)
const
7027 { md.brick_call(ib, version, 0); }
7029 void theta_method_dispatcher::asm_complex_tangent_terms
7030 (
const model &md,
size_type ib, model::complex_matlist &,
7031 std::vector<model::complex_veclist> &,
7032 std::vector<model::complex_veclist> &,
7033 build_version version)
const
7034 { md.brick_call(ib, version, 0); }
7036 theta_method_dispatcher::theta_method_dispatcher(
const std::string &THETA)
7037 : virtual_dispatcher(2) {
7038 param_names.push_back(THETA);
7042 (
model &md, dal::bit_vector ibricks,
const std::string &THETA) {
7043 pdispatcher pdispatch = std::make_shared<theta_method_dispatcher>(THETA);
7044 for (dal::bv_visitor i(ibricks); !i.finished(); ++i)
7049 (
model &md,
const std::string &U,
const std::string &V,
7050 const std::string &pdt,
const std::string &ptheta) {
7056 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for time step");
7058 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for parameter theta");
7061 scalar_type(1) - scalar_type(1) / theta[0]),
7064 scalar_type(1) / (theta[0]*dt[0])),
7067 -scalar_type(1) / (theta[0]*dt[0])),
7071 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for time step");
7072 const model_real_plain_vector &theta = md.
real_variable(ptheta);
7073 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for parameter theta");
7076 scalar_type(1) - scalar_type(1) / theta[0]),
7079 scalar_type(1) / (theta[0]*dt[0])),
7082 -scalar_type(1) / (theta[0]*dt[0])),
7094 (
model &md,
size_type id2dt2b,
const std::string &U,
const std::string &V,
7095 const std::string &pdt,
const std::string &ptwobeta,
7096 const std::string &pgamma) {
7106 if (twobeta != gamma) {
7108 md.set_dispatch_coeff();
7112 md.
assembly(model::BUILD_RHS_WITH_LIN);
7115 model_complex_plain_vector W(nbdof), RHS(nbdof);
7116 gmm::copy(gmm::sub_vector(md.
complex_rhs(), md.interval_of_variable(U)),
7121 gmm::cg(md.linear_complex_matrix_term(id2dt2b, 0),
7122 W, RHS, gmm::identity_matrix(), iter);
7123 GMM_ASSERT1(iter.converged(),
"Velocity not well computed");
7125 gmm::scaled(W, complex_type(1)/(twobeta*dt)),
7129 if (twobeta != gamma) {
7131 md.set_dispatch_coeff();
7135 GMM_ASSERT1(
false,
"to be done");
7144 if (twobeta != gamma) {
7146 md.set_dispatch_coeff();
7150 md.
assembly(model::BUILD_RHS_WITH_LIN);
7153 model_real_plain_vector W(nbdof), RHS(nbdof);
7154 gmm::copy(gmm::sub_vector(md.
real_rhs(), md.interval_of_variable(U)),
7159 gmm::cg(md.linear_real_matrix_term(id2dt2b, 0),
7160 W, RHS, gmm::identity_matrix(), iter);
7161 GMM_ASSERT1(iter.converged(),
"Velocity not well computed");
7163 gmm::scaled(W, scalar_type(1)/(twobeta*dt)),
7167 if (twobeta != gamma) {
7169 md.set_dispatch_coeff();
7184 class midpoint_dispatcher :
public virtual_dispatcher {
7186 gmm::uint64_type id_num;
7190 typedef model::build_version build_version;
7192 void set_dispatch_coeff(
const model &md,
size_type ib)
const {
7193 md.matrix_coeff_of_brick(ib) = scalar_type(1)/scalar_type(2);
7194 md.rhs_coeffs_of_brick(ib)[0] = scalar_type(1);
7195 md.rhs_coeffs_of_brick(ib)[1] = scalar_type(1)/scalar_type(2);
7198 template <
typename MATLIST,
typename VECTLIST>
7199 inline void next_iter(
const model &md,
size_type ib,
7200 const model::varnamelist &vl,
7201 const model::varnamelist &dl,
7203 VECTLIST &vectl, VECTLIST &vectl_sym,
7204 bool first_iter)
const {
7206 pbrick pbr = md.brick_pointer(ib);
7210 if (!(pbr->is_linear()))
7211 md.add_temporaries(vl, id_num);
7212 md.add_temporaries(dl, id_num);
7217 if (pbr->is_linear()) {
7220 if (first_iter) md.update_brick(ib, model::BUILD_RHS);
7223 md.linear_brick_add_to_rhs(ib, 1, 0);
7228 (
const model &md,
size_type ib,
const model::varnamelist &vl,
7229 const model::varnamelist &dl, model::real_matlist &matl,
7230 std::vector<model::real_veclist> &vectl,
7231 std::vector<model::real_veclist> &vectl_sym,
bool first_iter)
const {
7232 next_iter(md, ib, vl, dl, matl, vectl, vectl_sym, first_iter);
7235 void next_complex_iter
7236 (
const model &md,
size_type ib,
const model::varnamelist &vl,
7237 const model::varnamelist &dl,
7238 model::complex_matlist &matl,
7239 std::vector<model::complex_veclist> &vectl,
7240 std::vector<model::complex_veclist> &vectl_sym,
7241 bool first_iter)
const {
7242 next_iter(md, ib, vl, dl, matl, vectl, vectl_sym, first_iter);
7245 void asm_real_tangent_terms
7246 (
const model &md,
size_type ib, model::real_matlist &,
7247 std::vector<model::real_veclist> &vectl,
7248 std::vector<model::real_veclist> &vectl_sym,
7249 build_version version)
const {
7251 scalar_type half = scalar_type(1)/scalar_type(2);
7252 pbrick pbr = md.brick_pointer(ib);
7255 const model::varnamelist &vl = md.varnamelist_of_brick(ib);
7256 const model::varnamelist &dl = md.datanamelist_of_brick(ib);
7258 if (!(pbr->is_linear())) {
7259 for (
size_type i = 0; i < vl.size(); ++i) {
7260 bool is_uptodate = md.temporary_uptodate(vl[i], id_num, ind);
7261 if (!is_uptodate && ind !=
size_type(-1))
7262 gmm::add(gmm::scaled(md.real_variable(vl[i], 0), half),
7263 gmm::scaled(md.real_variable(vl[i], 1), half),
7264 md.set_real_variable(vl[i], ind));
7265 md.set_default_iter_of_variable(vl[i], ind);
7270 for (
size_type i = 0; i < dl.size(); ++i) {
7271 bool is_uptodate = md.temporary_uptodate(dl[i], id_num, ind);
7272 if (!is_uptodate && ind !=
size_type(-1)) {
7273 gmm::add(gmm::scaled(md.real_variable(dl[i], 0), half),
7274 gmm::scaled(md.real_variable(dl[i], 1), half),
7275 md.set_real_variable(dl[i], ind));
7277 md.set_default_iter_of_variable(dl[i], ind);
7281 md.brick_call(ib, version, 0);
7282 if (pbr->is_linear()) {
7286 md.linear_brick_add_to_rhs(ib, 1, 1);
7289 md.reset_default_iter_of_variables(dl);
7290 if (!(pbr->is_linear()))
7291 md.reset_default_iter_of_variables(vl);
7294 virtual void asm_complex_tangent_terms
7295 (
const model &md,
size_type ib, model::complex_matlist &,
7296 std::vector<model::complex_veclist> &vectl,
7297 std::vector<model::complex_veclist> &vectl_sym,
7298 build_version version)
const {
7300 scalar_type half = scalar_type(1)/scalar_type(2);
7301 pbrick pbr = md.brick_pointer(ib);
7304 const model::varnamelist &vl = md.varnamelist_of_brick(ib);
7305 const model::varnamelist &dl = md.datanamelist_of_brick(ib);
7307 if (!(pbr->is_linear())) {
7308 for (
size_type i = 0; i < vl.size(); ++i) {
7309 bool is_uptodate = md.temporary_uptodate(vl[i], id_num, ind);
7310 if (!is_uptodate && ind !=
size_type(-1))
7311 gmm::add(gmm::scaled(md.complex_variable(vl[i], 0), half),
7312 gmm::scaled(md.complex_variable(vl[i], 1), half),
7313 md.set_complex_variable(vl[i], ind));
7314 md.set_default_iter_of_variable(vl[i], ind);
7319 for (
size_type i = 0; i < dl.size(); ++i) {
7320 bool is_uptodate = md.temporary_uptodate(dl[i], id_num, ind);
7321 if (!is_uptodate && ind !=
size_type(-1)) {
7322 gmm::add(gmm::scaled(md.complex_variable(dl[i], 0), half),
7323 gmm::scaled(md.complex_variable(dl[i], 1), half),
7324 md.set_complex_variable(dl[i], ind));
7326 md.set_default_iter_of_variable(dl[i], ind);
7330 md.brick_call(ib, version, 0);
7331 if (pbr->is_linear()) {
7335 md.linear_brick_add_to_rhs(ib, 1, 1);
7338 md.reset_default_iter_of_variables(dl);
7339 if (!(pbr->is_linear()))
7340 md.reset_default_iter_of_variables(vl);
7343 midpoint_dispatcher() : virtual_dispatcher(2)
7344 { id_num = act_counter(); }
7349 pdispatcher pdispatch = std::make_shared<midpoint_dispatcher>();
7350 for (dal::bv_visitor i(ibricks); !i.finished(); ++i)
Takes a matrix or vector, or vector of matrices or vectors and creates an empty copy on each thread.
bool context_check() const
return true if update_from_context was called
im_data provides indexing to the integration points of a mesh im object.
Describe a finite element method linked to a mesh.
virtual dim_type get_qdim() const
Return the Q dimension.
virtual size_type nb_dof() const
Return the total number of degrees of freedom.
const mesh & linked_mesh() const
Return a reference to the underlying mesh.
Describe an integration method linked to a mesh.
const mesh & linked_mesh() const
Give a reference to the linked mesh of type mesh.
structure used to hold a set of convexes and/or convex faces.
int region_is_faces_of(const getfem::mesh &m1, const mesh_region &rg2, const getfem::mesh &m2) const
Test if the region is a boundary of a list of faces of elements of region rg.
Describe a mesh (collection of convexes (elements) and points).
const mesh_region region(size_type id) const
Return the region of index 'id'.
`‘Model’' variables store the variables, the data and the description of a model.
void add_initialized_matrix_data(const std::string &name, const base_matrix &M)
Add a fixed size data (assumed to be a matrix) to the model and initialized with M.
size_type nb_dof(bool with_internal=false) const
Total number of degrees of freedom in the model.
void delete_brick(size_type ib)
Delete the brick of index ib from the model.
void add_affine_dependent_variable(const std::string &name, const std::string &org_name, scalar_type alpha=scalar_type(1))
Add a "virtual" variable be an affine depedent variable with respect to another variable.
void add_macro(const std::string &name, const std::string &expr)
Add a macro definition for the high generic assembly language.
void disable_variable(const std::string &name)
Disable a variable (and its attached mutlipliers).
bool is_true_data(const std::string &name) const
States if a name corresponds to a declared data.
void change_mims_of_brick(size_type ib, const mimlist &ml)
Change the mim list of a brick.
void add_filtered_fem_variable(const std::string &name, const mesh_fem &mf, size_type region, size_type niter=1)
Add a variable linked to a fem with the dof filtered with respect to a mesh region.
virtual void assembly(build_version version)
Assembly of the tangent system taking into account all enabled terms in the model.
void delete_variable(const std::string &varname)
Delete a variable or data of the model.
const std::string & varname_of_brick(size_type ind_brick, size_type ind_var)
Gives the name of the variable of index ind_var of the brick of index ind_brick.
void add_fixed_size_data(const std::string &name, size_type size, size_type niter=1)
Add a fixed size data to the model.
bool is_linear() const
Return true if all the model terms are linear.
void add_fem_variable(const std::string &name, const mesh_fem &mf, size_type niter=1)
Add a variable being the dofs of a finite element method to the model.
const model_complex_plain_vector & complex_variable(const std::string &name, size_type niter) const
Gives the access to the vector value of a variable.
bool is_data(const std::string &name) const
States if a name corresponds to a declared data or disabled variable.
size_type add_brick(pbrick pbr, const varnamelist &varnames, const varnamelist &datanames, const termlist &terms, const mimlist &mims, size_type region)
Add a brick to the model.
void enable_brick(size_type ib)
Enable a brick.
void listvar(std::ostream &ost) const
List the model variables and constant.
const mesh_fem & mesh_fem_of_variable(const std::string &name) const
Gives the access to the mesh_fem of a variable if any.
const std::string & dataname_of_brick(size_type ind_brick, size_type ind_data)
Gives the name of the data of index ind_data of the brick of index ind_brick.
const model_real_plain_vector & real_rhs(bool with_internal=false) const
Gives access to the right hand side of the tangent linear system.
void add_im_data(const std::string &name, const im_data &imd, size_type niter=1)
Add data defined at integration points.
void add_multiplier(const std::string &name, const mesh_fem &mf, const std::string &primal_name, size_type niter=1)
Add a particular variable linked to a fem being a multiplier with respect to a primal variable.
void change_data_of_brick(size_type ib, const varnamelist &vl)
Change the data list of a brick.
const mesh_fem * pmesh_fem_of_variable(const std::string &name) const
Gives a pointer to the mesh_fem of a variable if any.
void resize_fixed_size_variable(const std::string &name, size_type size)
Resize a fixed size variable (or data) of the model.
bool macro_exists(const std::string &name) const
Says if a macro of that name has been defined.
model_real_plain_vector & set_real_variable(const std::string &name, size_type niter) const
Gives the write access to the vector value of a variable.
const model_complex_plain_vector & complex_rhs() const
Gives access to the right hand side of the tangent linear system.
void add_mim_to_brick(size_type ib, const mesh_im &mim)
Add an integration method to a brick.
void disable_brick(size_type ib)
Disable a brick.
bool is_disabled_variable(const std::string &name) const
States if a variable is disabled (treated as data).
void change_terms_of_brick(size_type ib, const termlist &terms)
Change the term list of a brick.
const model_real_plain_vector & real_variable(const std::string &name, size_type niter) const
Gives the access to the vector value of a variable.
bool is_complex() const
Boolean which says if the model deals with real or complex unknowns and data.
std::string Neumann_term(const std::string &varname, size_type region)
Gives the assembly string corresponding to the Neumann term of the fem variable varname on region.
bool has_internal_variables() const
Return true if the model has at least one internal variable.
void listbricks(std::ostream &ost, size_type base_id=0) const
List the model bricks.
bool is_internal_variable(const std::string &name) const
States if a variable is condensed out of the global system.
void add_time_dispatcher(size_type ibrick, pdispatcher pdispatch)
Add a time dispacther to a brick.
void add_fixed_size_variable(const std::string &name, size_type size, size_type niter=1)
Add a fixed size variable to the model assumed to be a vector.
void add_im_variable(const std::string &name, const im_data &imd, size_type niter=1)
Add variable defined at integration points.
void add_interpolate_transformation(const std::string &name, pinterpolate_transformation ptrans)
Add an interpolate transformation to the model to be used with the generic assembly.
void check_brick_stiffness_rhs(size_type ind_brick) const
check consistency of RHS and Stiffness matrix for brick with
model_complex_plain_vector & set_complex_variable(const std::string &name, size_type niter) const
Gives the write access to the vector value of a variable.
void change_variables_of_brick(size_type ib, const varnamelist &vl)
Change the variable list of a brick.
virtual void next_iter()
For transient problems.
void enable_variable(const std::string &name, bool enabled=true)
Enable a variable (and its attached mutlipliers).
void add_fem_data(const std::string &name, const mesh_fem &mf, dim_type qdim=1, size_type niter=1)
Add a data being the dofs of a finite element method to the model.
void change_update_flag_of_brick(size_type ib, bool flag)
Change the update flag of a brick.
std::string new_name(const std::string &name)
Gives a non already existing variable name begining by name.
void touch_brick(size_type ib)
Force the re-computation of a brick for the next assembly.
virtual void first_iter()
For transient problems.
void add_initialized_tensor_data(const std::string &name, const base_tensor &t)
Add a fixed size data (assumed to be a tensor) to the model and initialized with t.
void add_internal_im_variable(const std::string &name, const im_data &imd)
Add internal variable, defined at integration points and condensed.
bool variable_exists(const std::string &name) const
States if a name corresponds to a declared variable.
void del_macro(const std::string &name)
Delete a previously defined macro definition.
The virtual brick has to be derived to describe real model bricks.
virtual void asm_real_tangent_terms(const model &, size_type, const model::varnamelist &, const model::varnamelist &, const model::mimlist &, model::real_matlist &, model::real_veclist &, model::real_veclist &, size_type, build_version) const
Assembly of bricks real tangent terms.
void check_stiffness_matrix_and_rhs(const model &, size_type, const model::termlist &tlist, const model::varnamelist &, const model::varnamelist &, const model::mimlist &, model::real_matlist &, model::real_veclist &, model::real_veclist &, size_type rg, const scalar_type delta=1e-8) const
check consistency of stiffness matrix and rhs
virtual void real_pre_assembly_in_serial(const model &, size_type, const model::varnamelist &, const model::varnamelist &, const model::mimlist &, model::real_matlist &, model::real_veclist &, model::real_veclist &, size_type, build_version) const
Peform any pre assembly action for real term assembly.
virtual void real_post_assembly_in_serial(const model &, size_type, const model::varnamelist &, const model::varnamelist &, const model::mimlist &, model::real_matlist &, model::real_veclist &, model::real_veclist &, size_type, build_version) const
Peform any post assembly action for real terms.
The Iteration object calculates whether the solution has reached the desired accuracy,...
Distribution of assembly results (matrices/vectors) for parallel assembly.
Miscelleanous assembly routines for common terms. Use the low-level generic assembly....
Compute the gradient of a field on a getfem::mesh_fem.
A language for generic assembly of pde boundary value problems.
Compilation and execution operations.
Interpolation of fields from a mesh_fem onto another.
Model representation in Getfem.
#define GETFEM_OMP_PARALLEL(body)
Organizes a proper parallel omp section:
void mult_add(const L1 &l1, const L2 &l2, L3 &l3)
*/
void copy(const L1 &l1, L2 &l2)
*/
void clear(L &l)
clear (fill with zeros) a vector or matrix.
void resize(V &v, size_type n)
*/
void mult(const L1 &l1, const L2 &l2, L3 &l3)
*/
void add(const L1 &l1, L2 &l2)
*/
computation of the condition number of dense matrices.
Extract a basis of the range of a (large sparse) matrix from the columns of this matrix.
Conjugate gradient iterative solver.
void asm_mass_matrix(const MAT &M, const mesh_im &mim, const mesh_fem &mf1, const mesh_region &rg=mesh_region::all_convexes())
generic mass matrix assembly (on the whole mesh or on the specified convex set or boundary)
void asm_stiffness_matrix_for_homogeneous_laplacian(const MAT &M, const mesh_im &mim, const mesh_fem &mf, const mesh_region &rg=mesh_region::all_convexes())
assembly of .
void asm_stiffness_matrix_for_laplacian(MAT &M, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data, const VECT &A, const mesh_region &rg=mesh_region::all_convexes())
assembly of , where is scalar.
void asm_stiffness_matrix_for_scalar_elliptic(MAT &M, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data, const VECT &A, const mesh_region &rg=mesh_region::all_convexes())
assembly of , where is a (symmetric positive definite) NxN matrix.
void asm_source_term(const VECT1 &B, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data, const VECT2 &F, const mesh_region &rg=mesh_region::all_convexes())
source term (for both volumic sources and boundary (Neumann) sources).
void asm_homogeneous_normal_source_term(VECT1 &B, const mesh_im &mim, const mesh_fem &mf, const VECT2 &F, const mesh_region &rg)
Homogeneous normal source term (for boundary (Neumann) condition).
void asm_homogeneous_Helmholtz(MAT &M, const mesh_im &mim, const mesh_fem &mf_u, const VECT &K_squared, const mesh_region &rg=mesh_region::all_convexes())
assembly of the term , for the helmholtz equation ( , with ).
void asm_lumped_mass_matrix_for_first_order_param(MAT &M, const mesh_im &mim, const mesh_fem &mf_u, const mesh_fem &mf_data, const VECT &F, const mesh_region &rg=mesh_region::all_convexes())
lumped mass matrix assembly with an additional parameter (on the whole mesh or on the specified bound...
void asm_Helmholtz(MAT &M, const mesh_im &mim, const mesh_fem &mf_u, const mesh_fem &mf_data, const VECT &K_squared, const mesh_region &rg=mesh_region::all_convexes())
assembly of the term , for the helmholtz equation ( , with ).
void asm_stiffness_matrix_for_homogeneous_laplacian_componentwise(const MAT &M, const mesh_im &mim, const mesh_fem &mf, const mesh_region &rg=mesh_region::all_convexes())
assembly of .
void asm_normal_source_term(VECT1 &B, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data, const VECT2 &F, const mesh_region &rg)
Normal source term (for boundary (Neumann) condition).
void asm_stokes_B(const MAT &B, const mesh_im &mim, const mesh_fem &mf_u, const mesh_fem &mf_p, const mesh_region &rg=mesh_region::all_convexes())
Build the mixed pressure term .
void asm_lumped_mass_matrix_for_first_order(const MAT &M, const mesh_im &mim, const mesh_fem &mf1, const mesh_region &rg=mesh_region::all_convexes())
lumped mass matrix assembly (on the whole mesh or on the specified boundary)
void asm_mass_matrix_param(const MAT &M, const mesh_im &mim, const mesh_fem &mf1, const mesh_fem &mf2, const mesh_fem &mf_data, const VECT &A, const mesh_region &rg=mesh_region::all_convexes())
generic mass matrix assembly with an additional parameter (on the whole mesh or on the specified boun...
void asm_homogeneous_source_term(const VECT1 &B, const mesh_im &mim, const mesh_fem &mf, const VECT2 &F, const mesh_region &rg=mesh_region::all_convexes())
source term (for both volumic sources and boundary (Neumann) sources).
void asm_qu_term(MAT &M, const mesh_im &mim, const mesh_fem &mf_u, const mesh_fem &mf_d, const VECT &Q, const mesh_region &rg)
assembly of
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
size_t size_type
used as the common size type in the library
size_type alpha(short_type n, short_type d)
Return the value of which is the number of monomials of a polynomial of variables and degree .
GEneric Tool for Finite Element Methods.
size_type APIDECL add_nonlinear_term(model &md, const mesh_im &mim, const std::string &expr, size_type region=size_type(-1), bool is_sym=false, bool is_coercive=false, const std::string &brickname="")
Add a nonlinear term given by the weak form language expression expr which will be assembled in regio...
size_type APIDECL add_pointwise_constraints_with_penalization(model &md, const std::string &varname, scalar_type penalisation_coeff, const std::string &dataname_pt, const std::string &dataname_unitv=std::string(), const std::string &dataname_val=std::string())
Add some pointwise constraints on the variable varname thanks to a penalization.
size_type APIDECL add_isotropic_linearized_elasticity_pstrain_brick(model &md, const mesh_im &mim, const std::string &varname, const std::string &data_E, const std::string &data_nu, size_type region)
Linear elasticity brick ( ).
bool is_old(const std::string &name)
Does the variable have Old_ prefix.
size_type APIDECL add_nonlinear_twodomain_term(model &md, const mesh_im &mim, const std::string &expr, size_type region, const std::string &secondary_domain, bool is_sym=false, bool is_coercive=false, const std::string &brickname="")
Adds a nonlinear term given by a weak form language expression like add_nonlinear_term function but f...
size_type APIDECL add_normal_source_term_brick(model &md, const mesh_im &mim, const std::string &varname, const std::string &dataexpr, size_type region)
Add a source term on the variable varname on a boundary region.
size_type APIDECL add_lumped_mass_for_first_order_brick(model &md, const mesh_im &mim, const std::string &varname, const std::string &dataexpr_rho=std::string(), size_type region=size_type(-1))
Lumped mass brick for first order.
size_type APIDECL add_Dirichlet_condition_with_multipliers(model &md, const mesh_im &mim, const std::string &varname, const std::string &multname, size_type region, const std::string &dataname=std::string())
Add a Dirichlet condition on the variable varname and the mesh region region.
void asm_stiffness_matrix_for_homogeneous_scalar_elliptic_componentwise(MAT &M, const mesh_im &mim, const mesh_fem &mf, const VECT &A, const mesh_region &rg=mesh_region::all_convexes())
The same but with a constant matrix.
void APIDECL change_penalization_coeff(model &md, size_type ind_brick, scalar_type penalisation_coeff)
Change the penalization coefficient of a Dirichlet condition with penalization brick.
size_type APIDECL add_isotropic_linearized_elasticity_brick(model &md, const mesh_im &mim, const std::string &varname, const std::string &dataname_lambda, const std::string &dataname_mu, size_type region=size_type(-1), const std::string &dataname_preconstraint=std::string())
Linear elasticity brick ( ).
size_type APIDECL add_pointwise_constraints_with_given_multipliers(model &md, const std::string &varname, const std::string &multname, const std::string &dataname_pt, const std::string &dataname_unitv=std::string(), const std::string &dataname_val=std::string())
Add some pointwise constraints on the variable varname using a given multiplier multname.
size_type APIDECL add_basic_d2_on_dt2_brick(model &md, const mesh_im &mim, const std::string &varnameU, const std::string &datanameV, const std::string &dataname_dt, const std::string &dataname_alpha, const std::string &dataname_rho=std::string(), size_type region=size_type(-1))
Basic d2/dt2 brick ( ).
pinterpolate_transformation interpolate_transformation_neighbor_instance()
Create a new instance of a transformation corresponding to the interpolation on the neighbor element.
const auto PREFIX_OLD
A prefix to refer to the previous version of a variable.
size_type APIDECL add_mass_brick(model &md, const mesh_im &mim, const std::string &varname, const std::string &dataexpr_rho=std::string(), size_type region=size_type(-1))
Mass brick ( ).
void interpolation_von_mises_or_tresca(const getfem::mesh_fem &mf_u, const getfem::mesh_fem &mf_vm, const VEC1 &U, VEC2 &VM, const getfem::mesh_fem &mf_lambda, const VEC3 &lambda, const getfem::mesh_fem &mf_mu, const VEC3 &mu, bool tresca)
Compute the Von-Mises stress of a field (valid for linearized elasticity in 2D and 3D)
void asm_stiffness_matrix_for_scalar_elliptic_componentwise(MAT &M, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data, const VECT &A, const mesh_region &rg=mesh_region::all_convexes())
The same but on each component of mf when mf has a qdim > 1.
size_type APIDECL add_Laplacian_brick(model &md, const mesh_im &mim, const std::string &varname, size_type region=size_type(-1))
Add a Laplacian term on the variable varname (in fact with a minus : :math:-\text{div}(\nabla u)).
size_type APIDECL add_source_term(model &md, const mesh_im &mim, const std::string &expr, size_type region=size_type(-1), const std::string &brickname=std::string(), const std::string &directvarname=std::string(), const std::string &directdataname=std::string(), bool return_if_nonlin=false)
Add a source term given by the assembly string expr which will be assembled in region region and with...
size_type APIDECL add_generic_elliptic_brick(model &md, const mesh_im &mim, const std::string &varname, const std::string &dataexpr, size_type region=size_type(-1))
Add an elliptic term on the variable varname.
size_type APIDECL add_Dirichlet_condition_with_Nitsche_method(model &md, const mesh_im &mim, const std::string &varname, const std::string &Neumannterm, const std::string &datagamma0, size_type region, scalar_type theta=scalar_type(0), const std::string &datag=std::string())
Add a Dirichlet condition on the variable varname and the mesh region region.
size_type APIDECL add_Helmholtz_brick(model &md, const mesh_im &mim, const std::string &varname, const std::string &dataexpr, size_type region=size_type(-1))
Add a Helmoltz brick to the model.
void asm_stiffness_matrix_for_vector_elliptic(MAT &M, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data, const VECT &A, const mesh_region &rg=mesh_region::all_convexes())
Assembly of , where is a NxNxQxQ (symmetric positive definite) tensor defined on mf_data.
void APIDECL velocity_update_for_order_two_theta_method(model &md, const std::string &U, const std::string &V, const std::string &pdt, const std::string &ptheta)
Function which udpate the velocity $v^{n+1}$ after the computation of the displacement $u^{n+1}$ and ...
void APIDECL velocity_update_for_Newmark_scheme(model &md, size_type id2dt2b, const std::string &U, const std::string &V, const std::string &pdt, const std::string &ptwobeta, const std::string &pgamma)
Function which udpate the velocity $v^{n+1}$ after the computation of the displacement $u^{n+1}$ and ...
size_type APIDECL add_generalized_Dirichlet_condition_with_multipliers(model &md, const mesh_im &mim, const std::string &varname, const std::string &multname, size_type region, const std::string &dataname, const std::string &Hname)
Add a generalized Dirichlet condition on the variable varname and the mesh region region.
size_type APIDECL add_pointwise_constraints_with_multipliers(model &md, const std::string &varname, const std::string &dataname_pt, const std::string &dataname_unitv=std::string(), const std::string &dataname_val=std::string())
Add some pointwise constraints on the variable varname using multiplier.
void interpolation(const mesh_fem &mf_source, const mesh_fem &mf_target, const VECTU &U, VECTV &V, int extrapolation=0, double EPS=1E-10, mesh_region rg_source=mesh_region::all_convexes(), mesh_region rg_target=mesh_region::all_convexes())
interpolation/extrapolation of (mf_source, U) on mf_target.
size_type APIDECL add_twodomain_source_term(model &md, const mesh_im &mim, const std::string &expr, size_type region, const std::string &secondary_domain, const std::string &brickname=std::string(), const std::string &directvarname=std::string(), const std::string &directdataname=std::string(), bool return_if_nonlin=false)
Adds a source term given by a weak form language expression like add_source_term function but for an ...
void APIDECL compute_isotropic_linearized_Von_Mises_pstress(model &md, const std::string &varname, const std::string &data_E, const std::string &data_nu, const mesh_fem &mf_vm, model_real_plain_vector &VM)
Compute the Von-Mises stress of a displacement field for isotropic linearized elasticity in 3D or in ...
void asm_stiffness_matrix_for_laplacian_componentwise(MAT &M, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data, const VECT &A, const mesh_region &rg=mesh_region::all_convexes())
The same as getfem::asm_stiffness_matrix_for_laplacian , but on each component of mf when mf has a qd...
size_type APIDECL add_linear_incompressibility(model &md, const mesh_im &mim, const std::string &varname, const std::string &multname_pressure, size_type region=size_type(-1), const std::string &dataexpr_penal_coeff=std::string())
Mixed linear incompressibility condition brick.
std::shared_ptr< const virtual_brick > pbrick
type of pointer on a brick
const mesh_fem & classical_mesh_fem(const mesh &mesh, dim_type degree, dim_type qdim=1, bool complete=false)
Gives the descriptor of a classical finite element method of degree K on mesh.
size_type APIDECL add_generalized_Dirichlet_condition_with_penalization(model &md, const mesh_im &mim, const std::string &varname, scalar_type penalization_coeff, size_type region, const std::string &dataname, const std::string &Hname, const mesh_fem *mf_mult=0)
Add a Dirichlet condition on the variable varname and the mesh region region.
size_type APIDECL add_linear_term(model &md, const mesh_im &mim, const std::string &expr, size_type region=size_type(-1), bool is_sym=false, bool is_coercive=false, const std::string &brickname="", bool return_if_nonlin=false)
Add a term given by the weak form language expression expr which will be assembled in region region a...
size_type APIDECL add_linear_twodomain_term(model &md, const mesh_im &mim, const std::string &expr, size_type region, const std::string &secondary_domain, bool is_sym=false, bool is_coercive=false, const std::string &brickname="", bool return_if_nonlin=false)
Adds a linear term given by a weak form language expression like add_linear_term function but for an ...
void APIDECL compute_isotropic_linearized_Von_Mises_pstrain(model &md, const std::string &varname, const std::string &data_E, const std::string &data_nu, const mesh_fem &mf_vm, model_real_plain_vector &VM)
Compute the Von-Mises stress of a displacement field for isotropic linearized elasticity in 3D or in ...
size_type APIDECL add_generalized_Dirichlet_condition_with_Nitsche_method(model &md, const mesh_im &mim, const std::string &varname, const std::string &Neumannterm, const std::string &datagamma0, size_type region, scalar_type theta, const std::string &datag, const std::string &dataH)
Add a Dirichlet condition on the variable varname and the mesh region region.
size_type APIDECL add_Dirichlet_condition_with_penalization(model &md, const mesh_im &mim, const std::string &varname, scalar_type penalization_coeff, size_type region, const std::string &dataname=std::string(), const mesh_fem *mf_mult=0)
Add a Dirichlet condition on the variable varname and the mesh region region.
void APIDECL add_midpoint_dispatcher(model &md, dal::bit_vector ibricks)
Add a midpoint time dispatcher to a list of bricks.
size_type APIDECL add_Fourier_Robin_brick(model &md, const mesh_im &mim, const std::string &varname, const std::string &dataexpr, size_type region)
Add a Fourier-Robin brick to the model.
void APIDECL add_theta_method_dispatcher(model &md, dal::bit_vector ibricks, const std::string &THETA)
Add a theta-method time dispatcher to a list of bricks.
size_type APIDECL add_normal_Dirichlet_condition_with_penalization(model &md, const mesh_im &mim, const std::string &varname, scalar_type penalization_coeff, size_type region, const std::string &dataname=std::string(), const mesh_fem *mf_mult=0)
Add a Dirichlet condition to the normal component of the vector (or tensor) valued variable varname a...
size_type APIDECL add_isotropic_linearized_elasticity_pstress_brick(model &md, const mesh_im &mim, const std::string &varname, const std::string &data_E, const std::string &data_nu, size_type region)
Linear elasticity brick ( ).
size_type APIDECL add_basic_d_on_dt_brick(model &md, const mesh_im &mim, const std::string &varname, const std::string &dataname_dt, const std::string &dataname_rho=std::string(), size_type region=size_type(-1))
Basic d/dt brick ( ).
size_type APIDECL add_source_term_brick(model &md, const mesh_im &mim, const std::string &varname, const std::string &dataexpr, size_type region=size_type(-1), const std::string &directdataname=std::string())
Add a source term on the variable varname.
void asm_stiffness_matrix_for_homogeneous_scalar_elliptic(MAT &M, const mesh_im &mim, const mesh_fem &mf, const VECT &A, const mesh_region &rg=mesh_region::all_convexes())
The same but with a constant matrix.
size_type APIDECL add_normal_Dirichlet_condition_with_Nitsche_method(model &md, const mesh_im &mim, const std::string &varname, const std::string &Neumannterm, const std::string &datagamma0, size_type region, scalar_type theta=scalar_type(0), const std::string &datag=std::string())
Add a Dirichlet condition on the normal component of the variable varname and the mesh region region.
size_type APIDECL add_normal_Dirichlet_condition_with_multipliers(model &md, const mesh_im &mim, const std::string &varname, const std::string &multname, size_type region, const std::string &dataname=std::string())
Add a Dirichlet condition to the normal component of the vector (or tensor) valued variable varname a...
const APIDECL std::string & mult_varname_Dirichlet(model &md, size_type ind_brick)
When ind_brick is the index of a Dirichlet brick with multiplier on the model md, the function return...
size_type APIDECL add_Dirichlet_condition_with_simplification(model &md, const std::string &varname, size_type region, const std::string &dataname=std::string())
Add a (simple) Dirichlet condition on the variable varname and the mesh region region.
std::string no_old_prefix_name(const std::string &name)
Strip the variable name from prefix Old_ if it has one.
void asm_stiffness_matrix_for_homogeneous_vector_elliptic(MAT &M, const mesh_im &mim, const mesh_fem &mf, const VECT &A, const mesh_region &rg=mesh_region::all_convexes())
Assembly of , where is a NxNxQxQ (symmetric positive definite) constant tensor.