36 model::model(
bool comp_version) {
37 init(); complex_version = comp_version;
38 is_linear_ = is_symmetric_ = is_coercive_ =
true;
40 time_integration = 0; init_step =
false; time_step = scalar_type(1);
47 pstring s1 = std::make_shared<std::string>(
"Hess_u");
48 tree1.add_name(s1->c_str(), 6, 0, s1);
49 tree1.root->name =
"u";
50 tree1.root->op_type = GA_NAME;
51 tree1.root->node_type = GA_NODE_MACRO_PARAM;
53 tree1.root->nbc2 = ga_parse_prefix_operator(*s1);
54 tree1.root->nbc3 = ga_parse_prefix_test(*s1);
55 ga_macro gam1(
"Hess", tree1, 1);
56 macro_dict.add_macro(gam1);
59 pstring s2 = std::make_shared<std::string>(
"Div_u");
60 tree2.add_name(s2->c_str(), 5, 0, s2);
61 tree2.root->name =
"u";
62 tree2.root->op_type = GA_NAME;
63 tree2.root->node_type = GA_NODE_MACRO_PARAM;
65 tree2.root->nbc2 = ga_parse_prefix_operator(*s2);
66 tree2.root->nbc3 = ga_parse_prefix_test(*s2);
67 ga_macro gam2(
"Div", tree2, 1);
68 macro_dict.add_macro(gam2);
71 void model::var_description::set_size() {
73 v_num_var_iter.resize(n_iter);
74 v_num_iter.resize(n_iter);
75 size_type s = mf ? passociated_mf()->nb_dof()
76 : (imd ? imd->nb_filtered_index()
77 *imd->nb_tensor_elem()
82 complex_value[i].resize(s);
84 real_value[i].resize(s);
85 if (is_affine_dependent) {
87 affine_complex_value.resize(s);
89 affine_real_value.resize(s);
93 size_type model::var_description::add_temporary(gmm::uint64_type id_num) {
95 for (; nit < n_iter + n_temp_iter ; ++nit)
96 if (v_num_var_iter[nit] == id_num)
break;
97 if (nit >= n_iter + n_temp_iter) {
99 v_num_var_iter.resize(nit+1);
100 v_num_var_iter[nit] = id_num;
101 v_num_iter.resize(nit+1);
105 complex_value.resize(n_iter + n_temp_iter);
106 complex_value[nit].resize(s);
109 real_value.resize(n_iter + n_temp_iter);
110 real_value[nit].resize(s);
116 void model::var_description::clear_temporaries() {
120 complex_value.resize(n_iter);
122 real_value.resize(n_iter);
125 bool model::check_name_validity(
const std::string &name,
bool assert)
const {
127 if (variables.count(name) != 0) {
128 GMM_ASSERT1(!assert,
"Variable " << name <<
" already exists");
130 }
else if (variable_groups.count(name) != 0) {
132 name <<
" corresponds to an already existing group of "
137 name <<
" corresponds to an already existing macro");
139 }
else if (name.compare(
"X") == 0) {
140 GMM_ASSERT1(!assert,
"X is a reserved keyword of the generic "
141 "assembly language");
145 int ga_valid = ga_check_name_validity(name);
147 GMM_ASSERT1(!assert,
"Invalid variable name, corresponds to an "
148 "operator or function name of the generic assembly language");
150 }
else if (ga_valid == 2) {
151 GMM_ASSERT1(!assert,
"Invalid variable name having a reserved "
152 "prefix used by the generic assembly language");
154 }
else if (ga_valid == 3) {
155 std::string org_name = sup_previous_and_dot_to_varname(name);
156 if (org_name.size() < name.size() &&
157 variables.find(org_name) != variables.end()) {
159 "Dot and Previous are reserved prefix used for time "
160 "integration schemes");
165 bool valid = !name.empty() && isalpha(name[0]);
167 for (
size_type i = 1; i < name.size(); ++i)
168 if (!(isalnum(name[i]) || name[i] ==
'_')) valid =
false;
169 GMM_ASSERT1(!assert || valid,
170 "Illegal variable name : \"" << name <<
"\"");
175 std::string res_name = name;
176 bool valid = check_name_validity(res_name,
false);
177 for (
size_type i = 2; !valid && i < 50; ++i) {
179 m << name <<
'_' << i;
181 valid = check_name_validity(res_name,
false);
183 for (
size_type i = 2; !valid && i < 1000; ++i) {
185 m <<
"new_" << name <<
'_' << i;
187 valid = check_name_validity(res_name,
false);
189 GMM_ASSERT1(valid,
"Illegal variable name: " << name);
194 model::VAR_SET::const_iterator
195 model::find_variable(
const std::string &name)
const {
196 VAR_SET::const_iterator it = variables.find(name);
197 GMM_ASSERT1(it != variables.end(),
"Undefined variable " << name);
201 const model::var_description &
202 model::variable_description(
const std::string &name)
const {
203 return find_variable(name)->second;
206 std::string sup_previous_and_dot_to_varname(std::string v) {
207 if (!(v.compare(0, 8,
"Previous")) && (v[8] ==
'_' || v[9] ==
'_')) {
208 v = v.substr((v[8] ==
'_') ? 9 : 10);
210 if (!(v.compare(0, 3,
"Dot")) && (v[3] ==
'_' || v[4] ==
'_')) {
211 v = v.substr((v[3] ==
'_') ? 4 : 5);
218 return name.substr(0, PREFIX_OLD_LENGTH) ==
PREFIX_OLD;
222 return is_old(name) ? name.substr(PREFIX_OLD_LENGTH) : name;
226 if (
is_old(name))
return false;
227 VAR_SET::const_iterator it = find_variable(name);
228 if (!(it->second.is_variable))
return false;
229 if (it->second.is_affine_dependent)
230 it = variables.find(it->second.org_name);
231 return it->second.is_disabled;
235 if (
is_old(name))
return true;
236 VAR_SET::const_iterator it = find_variable(name);
237 if (it->second.is_affine_dependent)
238 it = variables.find(it->second.org_name);
239 return !(it->second.is_variable) || it->second.is_disabled;
243 return is_old(name) || !(variable_description(name).is_variable);
247 if (
is_old(name))
return false;
248 const auto &var_descr = variable_description(name);
249 return var_descr.is_internal && var_descr.is_enabled();
252 bool model::is_affine_dependent_variable(
const std::string &name)
const {
253 return !(
is_old(name)) && variable_description(name).is_affine_dependent;
257 model::org_variable(
const std::string &name)
const {
258 GMM_ASSERT1(is_affine_dependent_variable(name),
259 "For affine dependent variables only");
260 return variable_description(name).org_name;
264 model::factor_of_variable(
const std::string &name)
const {
265 return variable_description(name).alpha;
268 void model::set_factor_of_variable(
const std::string &name, scalar_type a) {
269 VAR_SET::iterator it = variables.find(name);
270 GMM_ASSERT1(it != variables.end(),
"Undefined variable " << name);
271 if (it->second.alpha != a) {
272 it->second.alpha = a;
273 for (
auto &v_num : it->second.v_num_data) v_num = act_counter();
277 bool model::is_im_data(
const std::string &name)
const {
282 model::pim_data_of_variable(
const std::string &name)
const {
286 const gmm::uint64_type &
287 model::version_number_of_data_variable(
const std::string &name,
289 VAR_SET::const_iterator it = find_variable(name);
290 if (niter ==
size_type(-1)) niter = it->second.default_iter;
291 return it->second.v_num_data[niter];
296 if (act_size_to_be_done) actualize_sizes();
298 return gmm::vect_size(crhs);
299 else if (with_internal && gmm::vect_size(full_rrhs))
300 return gmm::vect_size(full_rrhs);
302 return gmm::vect_size(rrhs);
305 void model::resize_global_system()
const {
308 for (
auto &&v : variables)
309 if (v.second.is_variable) {
310 if (v.second.is_disabled)
311 v.second.I = gmm::sub_interval(0,0);
312 else if (!v.second.is_affine_dependent && !v.second.is_internal) {
313 v.second.I = gmm::sub_interval(full_size, v.second.size());
314 full_size += v.second.size();
319 for (
auto &&v : variables)
320 if (v.second.is_internal && v.second.is_enabled()) {
321 v.second.I = gmm::sub_interval(full_size, v.second.size());
322 full_size += v.second.size();
325 for (
auto &&v : variables)
326 if (v.second.is_affine_dependent) {
327 v.second.I = variables.find(v.second.org_name)->second.I;
331 if (complex_version) {
340 if (full_size > primary_size) {
342 gmm::resize(internal_rTM, full_size-primary_size, primary_size);
351 for (dal::bv_visitor ib(valid_bricks); !ib.finished(); ++ib)
352 for (
const term_description &term : bricks[ib].tlist)
353 if (term.is_global) {
354 bricks[ib].terms_to_be_computed =
true;
359 void model::actualize_sizes()
const {
361 bool actualized =
false;
362 getfem::local_guard lock = locks_.get_lock();
363 if (actualized)
return;
365 act_size_to_be_done =
false;
367 std::map<std::string, std::vector<std::string> > multipliers;
368 std::set<std::string> tobedone;
384 for (
auto &&v : variables) {
385 const std::string &vname = v.first;
386 var_description &vdescr = v.second;
387 if (vdescr.mf && !vdescr.is_affine_dependent) {
388 if ((vdescr.filter & VDESCRFILTER_CTERM)
389 || (vdescr.filter & VDESCRFILTER_INFSUP)) {
390 VAR_SET::iterator vfilt = variables.find(vdescr.filter_var);
391 GMM_ASSERT1(vfilt != variables.end(),
"The primal variable of the"
392 " multiplier does not exist : " << vdescr.filter_var);
393 GMM_ASSERT1(vfilt->second.mf,
"The primal variable of the "
394 "multiplier should be a fem variable");
395 multipliers[vdescr.filter_var].push_back(vname);
396 if (vdescr.v_num < vdescr.mf->version_number() ||
397 vdescr.v_num < vfilt->second.mf->version_number()) {
398 tobedone.insert(vdescr.filter_var);
401 switch (vdescr.filter) {
402 case VDESCRFILTER_NO:
403 if (vdescr.v_num < vdescr.mf->version_number()) {
405 vdescr.v_num = act_counter();
408 case VDESCRFILTER_REGION:
409 if (vdescr.v_num < vdescr.mf->version_number()) {
411 dor = vdescr.mf->dof_on_region(vdescr.filter_region);
412 vdescr.partial_mf->adapt(dor);
414 vdescr.v_num = act_counter();
422 && vdescr.v_num < vdescr.imd->version_number()) {
424 vdescr.v_num = act_counter();
428 for (
auto &&v : variables) {
429 var_description &vdescr = v.second;
430 if (vdescr.mf && !(vdescr.is_affine_dependent) &&
431 ((vdescr.filter & VDESCRFILTER_CTERM)
432 || (vdescr.filter & VDESCRFILTER_INFSUP))) {
433 if (tobedone.count(vdescr.filter_var)) {
438 dal::bit_vector alldof; alldof.add(0, vdescr.mf->nb_dof());
439 vdescr.partial_mf->adapt(alldof);
441 vdescr.v_num = act_counter();
446 resize_global_system();
448 for (
const std::string &vname : tobedone) {
455 const std::vector<std::string> &mults = multipliers[vname];
456 const var_description &vdescr = variable_description(vname);
458 gmm::col_matrix< gmm::rsvector<scalar_type> > MGLOB;
459 if (mults.size() > 1) {
461 for (
const std::string &mult : mults)
462 s += variable_description(mult).mf->nb_dof();
466 std::set<size_type> glob_columns;
467 for (
const std::string &multname : mults) {
468 var_description &multdescr = variables.find(multname)->second;
474 gmm::col_matrix< gmm::rsvector<scalar_type> >
475 MM(vdescr.associated_mf().nb_dof(), multdescr.mf->nb_dof());
476 bool termadded =
false;
478 if (multdescr.filter & VDESCRFILTER_CTERM) {
480 for (dal::bv_visitor ib(valid_bricks); !ib.finished(); ++ib) {
481 const brick_description &brick = bricks[ib];
483 bool cplx =
is_complex() && brick.pbr->is_complex();
485 if (brick.tlist.size() == 0) {
486 bool varc =
false, multc =
false;
487 for (
const std::string &var : brick.vlist) {
488 if (multname.compare(var) == 0) multc =
true;
489 if (vname.compare(var) == 0) varc =
true;
492 GMM_ASSERT1(!cplx,
"Sorry, not taken into account");
493 generic_expressions.clear();
494 brick.terms_to_be_computed =
true;
495 update_brick(ib, BUILD_MATRIX);
496 if (generic_expressions.size()) {
497 GMM_TRACE2(
"Generic assembly for actualize sizes");
500 accumulated_distro<decltype(rTM)> distro_rTM(rTM);
502 ga_workspace workspace(*
this);
503 for (
const auto &ge : generic_expressions)
504 workspace.add_expression(ge.expr, ge.mim, ge.region,
505 2, ge.secondary_domain);
506 workspace.set_assembled_matrix(distro_rTM);
507 workspace.assembly(2);
510 gmm::add(gmm::sub_matrix(rTM, vdescr.I, multdescr.I), MM);
512 (gmm::sub_matrix(rTM, multdescr.I, vdescr.I)), MM);
518 for (
size_type j = 0; j < brick.tlist.size(); ++j) {
519 const term_description &term = brick.tlist[j];
521 if (term.is_matrix_term) {
522 if (term.is_global) {
523 bool varc =
false, multc =
false;
524 for (
const std::string &var : brick.vlist) {
525 if (multname.compare(var) == 0) multc =
true;
526 if (vname.compare(var) == 0) varc =
true;
529 GMM_ASSERT1(!cplx,
"Sorry, not taken into account");
530 generic_expressions.clear();
532 brick.terms_to_be_computed =
true;
533 update_brick(ib, BUILD_MATRIX);
536 gmm::add(gmm::sub_matrix(brick.rmatlist[j],
537 multdescr.I, vdescr.I),
539 gmm::add(gmm::transposed(gmm::sub_matrix
541 vdescr.I, multdescr.I)),
545 }
else if (multname.compare(term.var1) == 0 &&
546 vname.compare(term.var2) == 0) {
548 brick.terms_to_be_computed =
true;
549 update_brick(ib, BUILD_MATRIX);
553 gmm::add(gmm::transposed(gmm::real_part(brick.cmatlist[j])),
556 gmm::add(gmm::transposed(brick.rmatlist[j]), MM);
559 }
else if (multname.compare(term.var2) == 0 &&
560 vname.compare(term.var1) == 0) {
562 brick.terms_to_be_computed =
true;
563 update_brick(ib, BUILD_MATRIX);
567 gmm::add(gmm::real_part(brick.cmatlist[j]), MM);
577 GMM_WARNING1(
"No term found to filter multiplier " << multname
578 <<
". Variable is cancelled");
579 }
else if (multdescr.filter & VDESCRFILTER_INFSUP) {
580 mesh_region rg(multdescr.filter_region);
581 multdescr.filter_mim->linked_mesh().intersect_with_mpi_region(rg);
583 *(multdescr.mf), rg);
586 MPI_SUM_SPARSE_MATRIX(MM);
591 std::set<size_type> columns;
592 gmm::range_basis(MM, columns);
593 if (columns.size() == 0)
594 GMM_WARNING1(
"Empty basis found for multiplier " << multname);
596 if (mults.size() > 1) {
599 gmm::sub_interval(0, vdescr.associated_mf().nb_dof()),
600 gmm::sub_interval(s, multdescr.mf->nb_dof())));
602 glob_columns.insert(s + icol);
603 s += multdescr.mf->nb_dof();
605 dal::bit_vector kept;
608 if (multdescr.filter & VDESCRFILTER_REGION)
609 kept &= multdescr.mf->dof_on_region(multdescr.filter_region);
610 multdescr.partial_mf->adapt(kept);
611 multdescr.set_size();
612 multdescr.v_num = act_counter();
620 if (mults.size() > 1) {
621 gmm::range_basis(MGLOB, glob_columns, 1E-12, gmm::col_major(),
true);
628 for (
const std::string &multname : mults) {
629 var_description &multdescr = variables.find(multname)->second;
630 dal::bit_vector kept;
631 size_type nbdof = multdescr.mf->nb_dof();
632 for (
const size_type &icol : glob_columns)
633 if (icol >= s && icol < s + nbdof) kept.add(icol-s);
634 if (multdescr.filter & VDESCRFILTER_REGION)
635 kept &= multdescr.mf->dof_on_region(multdescr.filter_region);
636 multdescr.partial_mf->adapt(kept);
637 multdescr.set_size();
638 multdescr.v_num = act_counter();
639 s += multdescr.mf->nb_dof();
647 resize_global_system();
659 if (variables.size() == 0)
660 ost <<
"Model with no variable nor data" << endl;
662 ost <<
"List of model variables and data:" << endl;
663 for (
int vartype=0; vartype < 3; ++vartype)
664 for (
const auto &v : variables) {
665 const var_description &vdescr = v.second;
666 bool is_variable = vdescr.is_variable;
669 if (!is_variable || is_disabled)
continue;
670 }
else if (vartype == 1) {
671 if (!is_disabled)
continue;
672 }
else if (vartype == 2) {
673 if (is_variable)
continue;
675 ost << (is_variable ?
"Variable " :
"Data ");
676 ost << std::setw(30) << std::left << v.first;
677 ost << std::setw(2) << std::right << vdescr.n_iter;
678 ost << ((vdescr.n_iter == 1) ?
" copy " :
" copies ");
679 ost << (vdescr.mf ?
"fem dependant " :
"constant size ");
680 ost << std::setw(8) << std::right << vdescr.size();
682 ost << ((vdescr.size() > 1) ?
" doubles." :
" double.");
683 ost << (is_disabled ?
"\t (disabled)" :
"\t ");
684 if (vdescr.imd != 0) ost <<
"\t (is im_data)";
685 if (vdescr.is_affine_dependent) ost <<
"\t (is affine dependent)";
688 for (
const auto &vargroup : variable_groups) {
689 ost <<
"Variable group " << std::setw(30) << std::left
691 if (vargroup.second.size()) {
693 for (
const std::string &vname : vargroup.second) {
694 ost << (first ?
" " :
", ") << vname;
699 ost <<
" empty" << endl;
704 void model::listresiduals(std::ostream &ost)
const {
706 if (variables.size() == 0)
707 ost <<
"Model with no variable nor data" << endl;
710 for (
const auto &v : variables) {
711 if (v.second.is_variable) {
712 const model_real_plain_vector &rhs = v.second.is_internal
714 const gmm::sub_interval &II = interval_of_variable(v.first);
715 scalar_type res = gmm::vect_norm2(gmm::sub_vector(rhs, II));
716 if (!firstvar) cout <<
", ";
717 ost <<
"res_" << v.first <<
"= " << std::setw(11) << res;
727 bgeot::multi_index sizes(1);
733 const bgeot::multi_index &sizes,
735 check_name_validity(name);
736 variables.emplace(name, var_description(
true,
is_complex(), 0, 0, niter));
737 variables[name].qdims = sizes;
738 act_size_to_be_done =
true;
739 variables[name].set_size();
740 GMM_ASSERT1(variables[name].qdim(),
741 "Variables of null size are not allowed");
746 bgeot::multi_index sizes(1);
752 const bgeot::multi_index &sizes) {
753 GMM_ASSERT1(variables[name].mf == 0,
754 "Cannot explicitly resize a fem variable or data");
755 GMM_ASSERT1(variables[name].imd == 0,
756 "Cannot explicitly resize an im variable or data");
757 variables[name].qdims = sizes;
758 variables[name].set_size();
763 bgeot::multi_index sizes(1);
769 const bgeot::multi_index &sizes,
771 check_name_validity(name);
772 variables.emplace(name, var_description(
false,
is_complex(), 0, 0, niter));
773 variables[name].qdims = sizes;
774 GMM_ASSERT1(variables[name].qdim(),
"Data of null size are not allowed");
775 variables[name].set_size();
779 const base_matrix &M) {
782 GMM_ASSERT1(!(
is_complex()),
"Sorry, complex version to be done");
787 const base_complex_matrix &M) {
790 GMM_ASSERT1(!(
is_complex()),
"Sorry, complex version to be done");
795 const base_tensor &t) {
797 GMM_ASSERT1(!(
is_complex()),
"Sorry, complex version to be done");
802 const base_complex_tensor &t) {
804 GMM_ASSERT1(!(
is_complex()),
"Sorry, complex version to be done");
810 check_name_validity(name);
811 variables.emplace(name,
812 var_description(
true,
is_complex(), 0, &imd, niter));
813 variables[name].set_size();
815 act_size_to_be_done =
true;
821 variables[name].is_internal =
true;
826 check_name_validity(name);
827 variables.emplace(name,
828 var_description(
false,
is_complex(), 0, &imd, niter));
829 variables[name].set_size();
835 check_name_validity(name);
836 variables.emplace(name, var_description(
true,
is_complex(), &mf, 0, niter,
838 variables[name].set_size();
840 act_size_to_be_done =
true;
841 leading_dim = std::max(leading_dim, mf.
linked_mesh().dim());
847 check_name_validity(name);
848 variables.emplace(name, var_description(
true,
is_complex(), &mf, 0, niter,
849 VDESCRFILTER_REGION, region));
850 variables[name].set_size();
851 act_size_to_be_done =
true;
856 const std::string &org_name,
858 check_name_validity(name);
859 VAR_SET::const_iterator it = find_variable(org_name);
860 GMM_ASSERT1(it->second.is_variable && !(it->second.is_affine_dependent),
861 "The original variable should be a variable");
862 variables.emplace(name, variables[org_name]);
863 variables[name].is_affine_dependent =
true;
864 variables[name].org_name = org_name;
865 variables[name].alpha = alpha;
866 variables[name].set_size();
871 bgeot::multi_index sizes(1);
877 const bgeot::multi_index &sizes,
size_type niter) {
878 check_name_validity(name);
879 variables.emplace(name, var_description(
false,
is_complex(), &mf, 0, niter,
881 variables[name].qdims = sizes;
882 GMM_ASSERT1(variables[name].qdim(),
"Data of null size are not allowed");
883 variables[name].set_size();
888 const std::string &primal_name,
890 check_name_validity(name);
891 variables.emplace(name, var_description(
true,
is_complex(), &mf, 0, niter,
894 variables[name].set_size();
895 act_size_to_be_done =
true;
900 size_type region,
const std::string &primal_name,
902 check_name_validity(name);
903 variables.emplace(name, var_description(
true,
is_complex(), &mf, 0, niter,
904 VDESCRFILTER_REGION_CTERM, region,
906 variables[name].set_size();
907 act_size_to_be_done =
true;
912 const std::string &primal_name,
915 check_name_validity(name);
916 variables.emplace(name, var_description(
true,
is_complex(), &mf, 0, niter,
917 VDESCRFILTER_INFSUP, region,
919 variables[name].set_size();
920 act_size_to_be_done =
true;
930 VAR_SET::iterator it = variables.find(name);
931 GMM_ASSERT1(it != variables.end(),
"Undefined variable " << name);
932 it->second.is_disabled = !enabled;
933 for (
auto &&v : variables) {
934 if (((v.second.filter & VDESCRFILTER_INFSUP) ||
935 (v.second.filter & VDESCRFILTER_CTERM))
936 && name.compare(v.second.filter_var) == 0) {
937 v.second.is_disabled = !enabled;
939 if (v.second.is_variable && v.second.is_affine_dependent
940 && name.compare(v.second.org_name) == 0)
941 v.second.is_disabled = !enabled;
943 if (!act_size_to_be_done) resize_global_system();
951 check_name_validity(name.substr(0, name.find(
"(")));
952 macro_dict.add_macro(name, expr);
956 { macro_dict.del_macro(name); }
959 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
960 valid_bricks.del(ib);
961 active_bricks.del(ib);
963 for (
size_type i = 0; i < bricks[ib].mims.size(); ++i) {
964 const mesh_im *mim = bricks[ib].mims[i];
966 for (dal::bv_visitor ibb(valid_bricks); !ibb.finished(); ++ibb) {
967 for (
size_type j = 0; j < bricks[ibb].mims.size(); ++j)
968 if (bricks[ibb].mims[j] == mim) found =
true;
970 for (
const auto &v : variables) {
971 if (v.second.mf && (v.second.filter & VDESCRFILTER_INFSUP) &&
972 mim == v.second.filter_mim) found =
true;
974 if (!found) sup_dependency(*mim);
977 is_linear_ = is_symmetric_ = is_coercive_ =
true;
978 for (dal::bv_visitor ibb(valid_bricks); !ibb.finished(); ++ibb) {
979 is_linear_ = is_linear_ && bricks[ibb].pbr->is_linear();
980 is_symmetric_ = is_symmetric_ && bricks[ibb].pbr->is_symmetric();
981 is_coercive_ = is_coercive_ && bricks[ibb].pbr->is_coercive();
983 bricks[ib] = brick_description();
987 for (dal::bv_visitor ibb(valid_bricks); !ibb.finished(); ++ibb) {
988 for (
const auto &vname : bricks[ibb].vlist)
989 GMM_ASSERT1(varname.compare(vname),
990 "Cannot delete a variable which is still used by a brick");
991 for (
const auto &dname : bricks[ibb].dlist)
992 GMM_ASSERT1(varname.compare(dname),
993 "Cannot delete a data which is still used by a brick");
996 VAR_SET::const_iterator it = find_variable(varname);
1001 for(VAR_SET::iterator it2 = variables.begin();
1002 it2 != variables.end(); ++it2) {
1003 if (it != it2 && it2->second.mf && mf == it2->second.mf)
1006 if (!found) sup_dependency(*mf);
1008 if (it->second.filter & VDESCRFILTER_INFSUP) {
1009 const mesh_im *mim = it->second.filter_mim;
1011 for (dal::bv_visitor ibb(valid_bricks); !ibb.finished(); ++ibb) {
1012 for (
size_type j = 0; j < bricks[ibb].mims.size(); ++j)
1013 if (bricks[ibb].mims[j] == mim) found =
true;
1015 for (VAR_SET::iterator it2 = variables.begin();
1016 it2 != variables.end(); ++it2) {
1017 if (it != it2 && it2->second.mf &&
1018 (it2->second.filter & VDESCRFILTER_INFSUP) &&
1019 mim == it2->second.filter_mim) found =
true;
1021 if (!found) sup_dependency(*mim);
1025 if (it->second.imd != 0) sup_dependency(*(it->second.imd));
1027 variables.erase(varname);
1028 act_size_to_be_done =
true;
1032 const varnamelist &datanames,
1033 const termlist &terms,
1034 const mimlist &mims,
size_type region) {
1035 size_type ib = valid_bricks.first_false();
1037 for (
size_type i = 0; i < terms.size(); ++i)
1038 if (terms[i].is_global && terms[i].is_matrix_term && pbr->is_linear())
1039 GMM_ASSERT1(
false,
"Global linear matrix terms are not allowed");
1041 if (ib == bricks.size())
1042 bricks.push_back(brick_description(pbr, varnames, datanames, terms,
1045 bricks[ib] = brick_description(pbr, varnames, datanames, terms,
1047 active_bricks.add(ib);
1048 valid_bricks.add(ib);
1055 "Impossible to add a complex brick to a real model");
1057 bricks[ib].cmatlist.resize(terms.size());
1058 bricks[ib].cveclist[0].resize(terms.size());
1059 bricks[ib].cveclist_sym[0].resize(terms.size());
1061 bricks[ib].rmatlist.resize(terms.size());
1062 bricks[ib].rveclist[0].resize(terms.size());
1063 bricks[ib].rveclist_sym[0].resize(terms.size());
1065 is_linear_ = is_linear_ && pbr->is_linear();
1066 is_symmetric_ = is_symmetric_ && pbr->is_symmetric();
1067 is_coercive_ = is_coercive_ && pbr->is_coercive();
1069 for (
const auto &vname : varnames)
1070 GMM_ASSERT1(variables.count(vname),
1071 "Undefined model variable " << vname);
1073 for (
const auto &dname : datanames)
1074 GMM_ASSERT1(variables.count(dname),
1075 "Undefined model data or variable " << dname);
1081 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
1083 bricks[ib].mims.push_back(&mim);
1084 add_dependency(mim);
1088 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
1090 bricks[ib].tlist = terms;
1091 if (
is_complex() && bricks[ib].pbr->is_complex()) {
1092 bricks.back().cmatlist.resize(terms.size());
1093 bricks.back().cveclist[0].resize(terms.size());
1094 bricks.back().cveclist_sym[0].resize(terms.size());
1096 bricks.back().rmatlist.resize(terms.size());
1097 bricks.back().rveclist[0].resize(terms.size());
1098 bricks.back().rveclist_sym[0].resize(terms.size());
1103 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
1105 bricks[ib].vlist = vl;
1106 for (
const auto &v : vl)
1107 GMM_ASSERT1(variables.count(v),
"Undefined model variable " << v);
1111 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
1113 bricks[ib].dlist = dl;
1114 for (
const auto &v : dl)
1115 GMM_ASSERT1(variables.count(v),
"Undefined model variable " << v);
1119 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
1121 bricks[ib].mims = ml;
1122 for (
const auto &mim : ml) add_dependency(*mim);
1126 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
1128 bricks[ib].is_update_brick = flag;
1131 void model::set_time(scalar_type t,
bool to_init) {
1132 static const std::string varname(
"t");
1133 VAR_SET::iterator it = variables.find(varname);
1134 if (it == variables.end()) {
1137 GMM_ASSERT1(it->second.size() == 1,
"Time data should be of size 1");
1139 if (it == variables.end() || to_init) {
1147 scalar_type model::get_time() {
1148 static const std::string varname(
"t");
1149 set_time(scalar_type(0),
false);
1156 void model::call_init_affine_dependent_variables(
int version) {
1157 for (VAR_SET::iterator it = variables.begin();
1158 it != variables.end(); ++it) {
1159 var_description &vdescr = it->second;
1160 if (vdescr.is_variable && vdescr.ptsc) {
1162 vdescr.ptsc->init_affine_dependent_variables_precomputation(*
this);
1164 vdescr.ptsc->init_affine_dependent_variables(*
this);
1169 void model::shift_variables_for_time_integration() {
1170 for (VAR_SET::iterator it = variables.begin();
1171 it != variables.end(); ++it)
1172 if (it->second.is_variable && it->second.ptsc)
1173 it->second.ptsc->shift_variables(*
this);
1176 void model::add_time_integration_scheme(
const std::string &varname,
1177 ptime_scheme ptsc) {
1178 VAR_SET::iterator it = variables.find(varname);
1179 GMM_ASSERT1(it != variables.end(),
"Undefined variable " << varname);
1180 GMM_ASSERT1(it->second.is_variable && !(it->second.is_affine_dependent),
1181 "Cannot apply an integration scheme to " << varname);
1182 it->second.ptsc = ptsc;
1189 time_integration = 1;
1192 void model::copy_init_time_derivative() {
1194 for (VAR_SET::iterator it = variables.begin();
1195 it != variables.end(); ++it)
1196 if (it->second.is_variable && it->second.ptsc) {
1198 std::string name_v, name_previous_v;
1199 it->second.ptsc->time_derivative_to_be_initialized(name_v,
1202 if (name_v.size()) {
1220 class APIDECL first_order_theta_method_scheme
1221 :
public virtual_time_scheme {
1223 std::string U, U0, V, V0;
1228 virtual void init_affine_dependent_variables(model &md)
const {
1229 scalar_type dt = md.get_time_step();
1230 scalar_type a = scalar_type(1)/(theta*dt);
1231 scalar_type b = (scalar_type(1)-theta)/theta;
1232 md.set_factor_of_variable(V, a);
1233 if (md.is_complex()) {
1234 gmm::add(gmm::scaled(md.complex_variable(U0), -complex_type(a)),
1235 gmm::scaled(md.complex_variable(V0), -complex_type(b)),
1236 md.set_complex_constant_part(V));
1239 gmm::add(gmm::scaled(md.real_variable(U0), -a),
1240 gmm::scaled(md.real_variable(V0), -b),
1241 md.set_real_constant_part(V));
1246 virtual void init_affine_dependent_variables_precomputation(model &md)
1248 scalar_type dt = md.get_time_step();
1249 md.set_factor_of_variable(V, scalar_type(1)/dt);
1250 if (md.is_complex()) {
1251 gmm::copy(gmm::scaled(md.complex_variable(U0), -complex_type(1)/dt),
1252 md.set_complex_constant_part(V));
1255 gmm::copy(gmm::scaled(md.real_variable(U0), -scalar_type(1)/dt),
1256 md.set_real_constant_part(V));
1260 virtual void time_derivative_to_be_initialized
1261 (std::string &name_v, std::string &name_previous_v)
const
1262 {
if (theta != scalar_type(1)) { name_v = V; name_previous_v = V0; } }
1264 virtual void shift_variables(model &md)
const {
1265 if (md.is_complex()) {
1266 gmm::copy(md.complex_variable(U), md.set_complex_variable(U0));
1267 gmm::copy(md.complex_variable(V), md.set_complex_variable(V0));
1269 gmm::copy(md.real_variable(U), md.set_real_variable(U0));
1270 gmm::copy(md.real_variable(V), md.set_real_variable(V0));
1275 first_order_theta_method_scheme(model &md, std::string varname,
1278 U0 =
"Previous_" + U;
1280 V0 =
"Previous_Dot_" + U;
1282 GMM_ASSERT1(theta > scalar_type(0) && theta <= scalar_type(1),
1283 "Invalid value of theta parameter for the theta-method");
1285 if (!(md.variable_exists(V)))
1286 md.add_affine_dependent_variable(V, U);
1287 const mesh_fem *mf = md.pmesh_fem_of_variable(U);
1288 size_type s = md.is_complex() ? gmm::vect_size(md.complex_variable(U))
1289 : gmm::vect_size(md.real_variable(U));
1292 if (!(md.variable_exists(U0))) md.add_fem_data(U0, *mf);
1293 if (!(md.variable_exists(V0))) md.add_fem_data(V0, *mf);
1295 if (!(md.variable_exists(U0))) md.add_fixed_size_data(U0, s);
1296 if (!(md.variable_exists(V0))) md.add_fixed_size_data(V0, s);
1303 void add_theta_method_for_first_order(model &md,
const std::string &varname,
1304 scalar_type theta) {
1306 = std::make_shared<first_order_theta_method_scheme>(md, varname,theta);
1307 md.add_time_integration_scheme(varname, ptsc);
1316 class APIDECL second_order_theta_method_scheme
1317 :
public virtual_time_scheme {
1319 std::string U, U0, V, V0, A, A0;
1325 virtual void init_affine_dependent_variables(model &md)
const {
1326 scalar_type dt = md.get_time_step();
1327 md.set_factor_of_variable(V, scalar_type(1)/(theta*dt));
1328 md.set_factor_of_variable(A, scalar_type(1)/(theta*theta*dt*dt));
1329 if (md.is_complex()) {
1330 gmm::add(gmm::scaled(md.complex_variable(U0),
1331 -complex_type(1)/(theta*dt)),
1332 gmm::scaled(md.complex_variable(V0),
1333 -(complex_type(1)-complex_type(theta))/theta),
1334 md.set_complex_constant_part(V));
1335 gmm::add(gmm::scaled(md.complex_variable(U0),
1336 -complex_type(1)/(theta*theta*dt*dt)),
1337 gmm::scaled(md.complex_variable(A0),
1338 -(complex_type(1)-complex_type(theta))/theta),
1339 md.set_complex_constant_part(A));
1340 gmm::add(gmm::scaled(md.complex_variable(V0),
1341 -complex_type(1)/(theta*theta*dt)),
1342 md.set_complex_constant_part(A));
1346 gmm::add(gmm::scaled(md.real_variable(U0),
1347 -scalar_type(1)/(theta*dt)),
1348 gmm::scaled(md.real_variable(V0),
1349 -(scalar_type(1)-theta)/theta),
1350 md.set_real_constant_part(V));
1351 gmm::add(gmm::scaled(md.real_variable(U0),
1352 -scalar_type(1)/(theta*theta*dt*dt)),
1353 gmm::scaled(md.real_variable(A0),
1354 -(scalar_type(1)-theta)/theta),
1355 md.set_real_constant_part(A));
1356 gmm::add(gmm::scaled(md.real_variable(V0),
1357 -scalar_type(1)/(theta*theta*dt)),
1358 md.set_real_constant_part(A));
1365 virtual void init_affine_dependent_variables_precomputation(model &md)
1367 scalar_type dt = md.get_time_step();
1368 md.set_factor_of_variable(V, scalar_type(1)/dt);
1369 md.set_factor_of_variable(A, scalar_type(1)/(dt*dt));
1370 if (md.is_complex()) {
1371 gmm::copy(gmm::scaled(md.complex_variable(U0),
1372 -complex_type(1)/dt),
1373 md.set_complex_constant_part(V));
1374 gmm::add(gmm::scaled(md.complex_variable(U0),
1375 -complex_type(1)/(dt*dt)),
1376 gmm::scaled(md.complex_variable(V0),
1377 -complex_type(1)/dt),
1378 md.set_complex_constant_part(A));
1380 gmm::copy(gmm::scaled(md.real_variable(U0),
1381 -scalar_type(1)/dt),
1382 md.set_real_constant_part(V));
1383 gmm::add(gmm::scaled(md.real_variable(U0),
1384 -scalar_type(1)/(dt*dt)),
1385 gmm::scaled(md.real_variable(V0),
1386 -scalar_type(1)/dt),
1387 md.set_real_constant_part(A));
1391 virtual void time_derivative_to_be_initialized
1392 (std::string &name_v, std::string &name_previous_v)
const
1393 {
if (theta != scalar_type(1)) { name_v = A; name_previous_v = A0; } }
1395 virtual void shift_variables(model &md)
const {
1396 if (md.is_complex()) {
1397 gmm::copy(md.complex_variable(U), md.set_complex_variable(U0));
1398 gmm::copy(md.complex_variable(V), md.set_complex_variable(V0));
1399 gmm::copy(md.complex_variable(A), md.set_complex_variable(A0));
1401 gmm::copy(md.real_variable(U), md.set_real_variable(U0));
1402 gmm::copy(md.real_variable(V), md.set_real_variable(V0));
1403 gmm::copy(md.real_variable(A), md.set_real_variable(A0));
1408 second_order_theta_method_scheme(model &md, std::string varname,
1411 U0 =
"Previous_" + U;
1413 V0 =
"Previous_Dot_" + U;
1415 A0 =
"Previous_Dot2_" + U;
1417 GMM_ASSERT1(theta > scalar_type(0) && theta <= scalar_type(1),
1418 "Invalid value of theta parameter for the theta-method");
1420 if (!(md.variable_exists(V)))
1421 md.add_affine_dependent_variable(V, U);
1422 if (!(md.variable_exists(A)))
1423 md.add_affine_dependent_variable(A, U);
1425 const mesh_fem *mf = md.pmesh_fem_of_variable(U);
1426 size_type s = md.is_complex() ? gmm::vect_size(md.complex_variable(U))
1427 : gmm::vect_size(md.real_variable(U));
1430 if (!(md.variable_exists(U0))) md.add_fem_data(U0, *mf);
1431 if (!(md.variable_exists(V0))) md.add_fem_data(V0, *mf);
1432 if (!(md.variable_exists(A0))) md.add_fem_data(A0, *mf);
1434 if (!(md.variable_exists(U0))) md.add_fixed_size_data(U0, s);
1435 if (!(md.variable_exists(V0))) md.add_fixed_size_data(V0, s);
1436 if (!(md.variable_exists(A0))) md.add_fixed_size_data(A0, s);
1443 void add_theta_method_for_second_order(model &md,
const std::string &varname,
1444 scalar_type theta) {
1445 ptime_scheme ptsc = std::make_shared<second_order_theta_method_scheme>
1447 md.add_time_integration_scheme(varname, ptsc);
1457 class APIDECL Newmark_scheme
1458 :
public virtual_time_scheme {
1460 std::string U, U0, V, V0, A, A0;
1461 scalar_type beta, gamma;
1466 virtual void init_affine_dependent_variables(model &md)
const {
1467 scalar_type dt = md.get_time_step();
1468 scalar_type a0 = scalar_type(1)/(beta*dt*dt), a1 = dt*a0;
1469 scalar_type a2 = (scalar_type(1) - scalar_type(2)*beta)
1470 / (scalar_type(2)*beta);
1471 scalar_type b0 = gamma/(beta*dt), b1 = (beta-gamma)/beta;
1472 scalar_type b2 = dt*(scalar_type(1)-gamma/(scalar_type(2)*beta));
1474 md.set_factor_of_variable(V, b0);
1475 md.set_factor_of_variable(A, a0);
1476 if (md.is_complex()) {
1477 gmm::add(gmm::scaled(md.complex_variable(U0), -complex_type(b0)),
1478 gmm::scaled(md.complex_variable(V0), complex_type(b1)),
1479 md.set_complex_constant_part(V));
1480 gmm::add(gmm::scaled(md.complex_variable(A0), complex_type(b2)),
1481 md.set_complex_constant_part(V));
1482 gmm::add(gmm::scaled(md.complex_variable(U0), -complex_type(a0)),
1483 gmm::scaled(md.complex_variable(V0), -complex_type(a1)),
1484 md.set_complex_constant_part(A));
1485 gmm::add(gmm::scaled(md.complex_variable(A0), -complex_type(a2)),
1486 md.set_complex_constant_part(A));
1488 gmm::add(gmm::scaled(md.real_variable(U0), -b0),
1489 gmm::scaled(md.real_variable(V0), b1),
1490 md.set_real_constant_part(V));
1491 gmm::add(gmm::scaled(md.real_variable(A0), b2),
1492 md.set_real_constant_part(V));
1493 gmm::add(gmm::scaled(md.real_variable(U0), -a0),
1494 gmm::scaled(md.real_variable(V0), -a1),
1495 md.set_real_constant_part(A));
1496 gmm::add(gmm::scaled(md.real_variable(A0), -a2),
1497 md.set_real_constant_part(A));
1504 virtual void init_affine_dependent_variables_precomputation(model &md)
1506 scalar_type dt = md.get_time_step();
1507 md.set_factor_of_variable(V, scalar_type(1)/dt);
1508 md.set_factor_of_variable(A, scalar_type(1)/(dt*dt));
1509 if (md.is_complex()) {
1510 gmm::copy(gmm::scaled(md.complex_variable(U0),
1511 -complex_type(1)/dt),
1512 md.set_complex_constant_part(V));
1513 gmm::add(gmm::scaled(md.complex_variable(U0),
1514 -complex_type(1)/(dt*dt)),
1515 gmm::scaled(md.complex_variable(V0),
1516 -complex_type(1)/dt),
1517 md.set_complex_constant_part(A));
1519 gmm::copy(gmm::scaled(md.real_variable(U0),
1520 -scalar_type(1)/dt),
1521 md.set_real_constant_part(V));
1522 gmm::add(gmm::scaled(md.real_variable(U0),
1523 -scalar_type(1)/(dt*dt)),
1524 gmm::scaled(md.real_variable(V0),
1525 -scalar_type(1)/dt),
1526 md.set_real_constant_part(A));
1530 virtual void time_derivative_to_be_initialized
1531 (std::string &name_v, std::string &name_previous_v)
const {
1532 if (beta != scalar_type(0.5) || gamma != scalar_type(1))
1533 { name_v = A; name_previous_v = A0; }
1536 virtual void shift_variables(model &md)
const {
1537 if (md.is_complex()) {
1538 gmm::copy(md.complex_variable(U), md.set_complex_variable(U0));
1539 gmm::copy(md.complex_variable(V), md.set_complex_variable(V0));
1540 gmm::copy(md.complex_variable(A), md.set_complex_variable(A0));
1542 gmm::copy(md.real_variable(U), md.set_real_variable(U0));
1543 gmm::copy(md.real_variable(V), md.set_real_variable(V0));
1544 gmm::copy(md.real_variable(A), md.set_real_variable(A0));
1549 Newmark_scheme(model &md, std::string varname,
1550 scalar_type be, scalar_type ga) {
1552 U0 =
"Previous_" + U;
1554 V0 =
"Previous_Dot_" + U;
1556 A0 =
"Previous_Dot2_" + U;
1557 beta = be; gamma = ga;
1558 GMM_ASSERT1(beta > scalar_type(0) && beta <= scalar_type(1)
1559 && gamma >= scalar_type(0.5) && gamma <= scalar_type(1),
1560 "Invalid parameter values for the Newmark scheme");
1562 if (!(md.variable_exists(V)))
1563 md.add_affine_dependent_variable(V, U);
1564 if (!(md.variable_exists(A)))
1565 md.add_affine_dependent_variable(A, U);
1567 const mesh_fem *mf = md.pmesh_fem_of_variable(U);
1568 size_type s = md.is_complex() ? gmm::vect_size(md.complex_variable(U))
1569 : gmm::vect_size(md.real_variable(U));
1572 if (!(md.variable_exists(U0))) md.add_fem_data(U0, *mf);
1573 if (!(md.variable_exists(V0))) md.add_fem_data(V0, *mf);
1574 if (!(md.variable_exists(A0))) md.add_fem_data(A0, *mf);
1576 if (!(md.variable_exists(U0))) md.add_fixed_size_data(U0, s);
1577 if (!(md.variable_exists(V0))) md.add_fixed_size_data(V0, s);
1578 if (!(md.variable_exists(A0))) md.add_fixed_size_data(A0, s);
1585 void add_Newmark_scheme(model &md,
const std::string &varname,
1586 scalar_type beta, scalar_type gamma) {
1587 ptime_scheme ptsc = std::make_shared<Newmark_scheme>
1588 (md, varname, beta, gamma);
1589 md.add_time_integration_scheme(varname, ptsc);
1598 class APIDECL Houbolt_scheme
1599 :
public virtual_time_scheme {
1601 std::string U, U01, U02, U03, V, A;
1606 virtual void init_affine_dependent_variables(model &md)
const {
1607 scalar_type dt = md.get_time_step();
1608 scalar_type a0 = scalar_type(2)/(dt*dt);
1609 scalar_type a1 = scalar_type(5)/(dt*dt);
1610 scalar_type a2 = scalar_type(4)/(dt*dt);
1611 scalar_type a3 = scalar_type(1)/(dt*dt);
1612 scalar_type b0 = scalar_type(11)/(scalar_type(6)*dt);
1613 scalar_type b1 = scalar_type(18)/(scalar_type(6)*dt);
1614 scalar_type b2 = scalar_type(9)/(scalar_type(6)*dt);
1615 scalar_type b3 = scalar_type(2)/(scalar_type(6)*dt);
1617 md.set_factor_of_variable(V, b0);
1618 md.set_factor_of_variable(A, a0);
1619 if (md.is_complex()) {
1620 gmm::add(gmm::scaled(md.complex_variable(U01), -complex_type(b1)),
1621 gmm::scaled(md.complex_variable(U02), complex_type(b2)),
1622 md.set_complex_constant_part(V));
1623 gmm::add(gmm::scaled(md.complex_variable(U03), -complex_type(b3)),
1624 md.set_complex_constant_part(V));
1625 gmm::add(gmm::scaled(md.complex_variable(U01), -complex_type(a1)),
1626 gmm::scaled(md.complex_variable(U02), complex_type(a2)),
1627 md.set_complex_constant_part(A));
1628 gmm::add(gmm::scaled(md.complex_variable(U03), -complex_type(a3)),
1629 md.set_complex_constant_part(A));
1631 gmm::add(gmm::scaled(md.real_variable(U01), -b1),
1632 gmm::scaled(md.real_variable(U02), b2),
1633 md.set_real_constant_part(V));
1634 gmm::add(gmm::scaled(md.real_variable(U03), -b3),
1635 md.set_real_constant_part(V));
1636 gmm::add(gmm::scaled(md.real_variable(U01), -a1),
1637 gmm::scaled(md.real_variable(U02), a2),
1638 md.set_real_constant_part(A));
1639 gmm::add(gmm::scaled(md.real_variable(U03), -a3),
1640 md.set_real_constant_part(A));
1644 virtual void init_affine_dependent_variables_precomputation(model &md)
1649 virtual void time_derivative_to_be_initialized
1650 (std::string &name_v, std::string &name_previous_v)
const {
1652 (void) name_previous_v;
1655 virtual void shift_variables(model &md)
const {
1656 if (md.is_complex()) {
1657 gmm::copy(md.complex_variable(U02), md.set_complex_variable(U03));
1658 gmm::copy(md.complex_variable(U01), md.set_complex_variable(U02));
1659 gmm::copy(md.complex_variable(U), md.set_complex_variable(U01));
1661 gmm::copy(md.real_variable(U02), md.set_real_variable(U03));
1662 gmm::copy(md.real_variable(U01), md.set_real_variable(U02));
1663 gmm::copy(md.real_variable(U), md.set_real_variable(U01));
1668 Houbolt_scheme(model &md, std::string varname) {
1670 U01 =
"Previous_" + U;
1671 U02 =
"Previous2_" + U;
1672 U03 =
"Previous3_" + U;
1676 if (!(md.variable_exists(V)))
1677 md.add_affine_dependent_variable(V, U);
1678 if (!(md.variable_exists(A)))
1679 md.add_affine_dependent_variable(A, U);
1681 const mesh_fem *mf = md.pmesh_fem_of_variable(U);
1682 size_type s = md.is_complex() ? gmm::vect_size(md.complex_variable(U))
1683 : gmm::vect_size(md.real_variable(U));
1686 if (!(md.variable_exists(U01))) md.add_fem_data(U01, *mf);
1687 if (!(md.variable_exists(U02))) md.add_fem_data(U02, *mf);
1688 if (!(md.variable_exists(U03))) md.add_fem_data(U03, *mf);
1690 if (!(md.variable_exists(U01))) md.add_fixed_size_data(U01, s);
1691 if (!(md.variable_exists(U02))) md.add_fixed_size_data(U02, s);
1692 if (!(md.variable_exists(U03))) md.add_fixed_size_data(U03, s);
1699 void add_Houbolt_scheme(model &md,
const std::string &varname) {
1700 ptime_scheme ptsc = std::make_shared<Houbolt_scheme>
1702 md.add_time_integration_scheme(varname, ptsc);
1714 GMM_ASSERT1(valid_bricks[ibrick],
"Inexistent brick");
1715 pbrick pbr = bricks[ibrick].pbr;
1717 bricks[ibrick].pdispatch = pdispatch;
1720 = std::max(
size_type(1), pdispatch->nbrhs());
1722 gmm::resize(bricks[ibrick].coeffs, nbrhs);
1725 bricks[ibrick].cveclist.resize(nbrhs);
1726 bricks[ibrick].cveclist_sym.resize(nbrhs);
1728 bricks[ibrick].cveclist[k] = bricks[ibrick].cveclist[0];
1729 bricks[ibrick].cveclist_sym[k] = bricks[ibrick].cveclist_sym[0];
1732 bricks[ibrick].rveclist.resize(nbrhs);
1733 bricks[ibrick].rveclist_sym.resize(nbrhs);
1735 bricks[ibrick].rveclist[k] = bricks[ibrick].rveclist[0];
1736 bricks[ibrick].rveclist_sym[k] = bricks[ibrick].rveclist_sym[0];
1743 GMM_ASSERT1(valid_bricks[ind_brick],
"Inexistent brick");
1744 GMM_ASSERT1(ind_var < bricks[ind_brick].vlist.size(),
1745 "Inexistent brick variable");
1746 return bricks[ind_brick].vlist[ind_var];
1751 GMM_ASSERT1(valid_bricks[ind_brick],
"Inexistent brick");
1752 GMM_ASSERT1(ind_data < bricks[ind_brick].dlist.size(),
1753 "Inexistent brick data");
1754 return bricks[ind_brick].dlist[ind_data];
1758 if (valid_bricks.card() == 0)
1759 ost <<
"Model with no bricks" << endl;
1761 ost <<
"List of model bricks:" << endl;
1762 for (dal::bv_visitor i(valid_bricks); !i.finished(); ++i) {
1763 ost <<
"Brick " << std::setw(3) << std::right << i + base_id
1764 <<
" " << std::setw(20) << std::right
1765 << bricks[i].pbr->brick_name();
1766 if (!(active_bricks[i])) ost <<
" (deactivated)";
1767 if (bricks[i].pdispatch) ost <<
" (dispatched)";
1768 ost << endl <<
" concerned variables: " << bricks[i].vlist[0];
1769 for (
size_type j = 1; j < bricks[i].vlist.size(); ++j)
1770 ost <<
", " << bricks[i].vlist[j];
1772 ost <<
" brick with " << bricks[i].tlist.size() <<
" term";
1773 if (bricks[i].tlist.size() > 1) ost <<
"s";
1782 void model::brick_init(
size_type ib, build_version version,
1784 const brick_description &brick = bricks[ib];
1785 bool cplx =
is_complex() && brick.pbr->is_complex();
1788 for (
size_type j = 0; j < brick.tlist.size(); ++j) {
1789 const term_description &term = brick.tlist[j];
1790 bool isg = term.is_global;
1792 gmm::vect_size(crhs) : gmm::vect_size(rrhs);
1793 size_type nbd1 = isg ? nbgdof : variables[term.var1].size();
1794 size_type nbd2 = isg ? nbgdof : (term.is_matrix_term ?
1795 variables[term.var2].size() : 0);
1796 if (term.is_matrix_term &&
1797 (brick.pbr->is_linear() || (version | BUILD_MATRIX))) {
1798 if (version | BUILD_ON_DATA_CHANGE) {
1800 gmm::resize(brick.cmatlist[j], nbd1, nbd2);
1802 gmm::resize(brick.rmatlist[j], nbd1, nbd2);
1805 brick.cmatlist[j] = model_complex_sparse_matrix(nbd1, nbd2);
1807 brick.rmatlist[j] = model_real_sparse_matrix(nbd1, nbd2);
1810 if (brick.pbr->is_linear() || (version | BUILD_RHS)) {
1811 for (
size_type k = 0; k < brick.nbrhs; ++k) {
1813 if (k == rhs_ind)
gmm::clear(brick.cveclist[k][j]);
1815 if (term.is_symmetric && term.var1.compare(term.var2)) {
1816 if (k == rhs_ind)
gmm::clear(brick.cveclist_sym[k][j]);
1820 if (k == rhs_ind)
gmm::clear(brick.rveclist[k][j]);
1822 if (term.is_symmetric && term.var1.compare(term.var2)) {
1823 if (k == rhs_ind)
gmm::clear(brick.rveclist_sym[k][j]);
1832 void model::post_to_variables_step(){}
1834 void model::brick_call(
size_type ib, build_version version,
1837 const brick_description &brick = bricks[ib];
1838 bool cplx =
is_complex() && brick.pbr->is_complex();
1840 brick_init(ib, version, rhs_ind);
1844 brick.pbr->complex_pre_assembly_in_serial(*
this, ib, brick.vlist,
1845 brick.dlist, brick.mims,
1847 brick.cveclist[rhs_ind],
1848 brick.cveclist_sym[rhs_ind],
1849 brick.region, version);
1854 accumulated_distro<complex_matlist> cmatlist(brick.cmatlist);
1855 accumulated_distro<complex_veclist> cveclist(brick.cveclist[rhs_ind]);
1856 accumulated_distro<complex_veclist> cveclist_sym(brick.cveclist_sym[rhs_ind]);
1860 brick.pbr->asm_complex_tangent_terms(*
this, ib, brick.vlist,
1861 brick.dlist, brick.mims,
1865 brick.region, version);
1868 brick.pbr->complex_post_assembly_in_serial(*
this, ib, brick.vlist,
1869 brick.dlist, brick.mims,
1871 brick.cveclist[rhs_ind],
1872 brick.cveclist_sym[rhs_ind],
1873 brick.region, version);
1875 if (brick.is_update_brick)
1877 for (
auto &&mat : brick.cmatlist)
1880 for (
auto &&vecs : brick.cveclist)
1881 for (
auto &&vec : vecs)
1884 for (
auto &&vecs : brick.cveclist_sym)
1885 for (
auto &&vec : vecs)
1891 brick.pbr->real_pre_assembly_in_serial(*
this, ib, brick.vlist,
1892 brick.dlist, brick.mims,
1894 brick.rveclist[rhs_ind],
1895 brick.rveclist_sym[rhs_ind],
1896 brick.region, version);
1899 accumulated_distro<real_matlist> rmatlist(brick.rmatlist);
1900 accumulated_distro<real_veclist> rveclist(brick.rveclist[rhs_ind]);
1901 accumulated_distro<real_veclist> rveclist_sym(brick.rveclist_sym[rhs_ind]);
1905 brick.pbr->asm_real_tangent_terms(*
this, ib, brick.vlist,
1906 brick.dlist, brick.mims,
1914 brick.pbr->real_post_assembly_in_serial(*
this, ib, brick.vlist,
1915 brick.dlist, brick.mims,
1917 brick.rveclist[rhs_ind],
1918 brick.rveclist_sym[rhs_ind],
1919 brick.region, version);
1921 if (brick.is_update_brick)
1923 for (
auto &&mat : brick.rmatlist)
1926 for (
auto &&vecs : brick.rveclist)
1927 for (
auto &&vec : vecs)
1930 for (
auto &&vecs : brick.rveclist_sym)
1931 for (
auto &&vec : vecs)
1938 void model::set_dispatch_coeff() {
1939 for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) {
1940 brick_description &brick = bricks[ib];
1941 if (brick.pdispatch)
1942 brick.pdispatch->set_dispatch_coeff(*
this, ib);
1948 context_check();
if (act_size_to_be_done) actualize_sizes();
1949 for (
auto && v : variables) v.second.clear_temporaries();
1951 set_dispatch_coeff();
1953 for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) {
1954 brick_description &brick = bricks[ib];
1955 if (brick.pdispatch) {
1957 brick.pdispatch->next_complex_iter(*
this, ib, brick.vlist,
1959 brick.cmatlist, brick.cveclist,
1960 brick.cveclist_sym,
true);
1962 brick.pdispatch->next_real_iter(*
this, ib, brick.vlist, brick.dlist,
1963 brick.rmatlist, brick.rveclist,
1964 brick.rveclist_sym,
true);
1970 context_check();
if (act_size_to_be_done) actualize_sizes();
1971 set_dispatch_coeff();
1973 for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) {
1974 brick_description &brick = bricks[ib];
1975 if (brick.pdispatch) {
1977 brick.pdispatch->next_complex_iter(*
this, ib, brick.vlist,
1979 brick.cmatlist, brick.cveclist,
1980 brick.cveclist_sym,
false);
1982 brick.pdispatch->next_real_iter(*
this, ib, brick.vlist, brick.dlist,
1983 brick.rmatlist, brick.rveclist,
1984 brick.rveclist_sym,
false);
1988 for (
auto &&v : variables)
1989 for (
size_type i = 1; i < v.second.n_iter; ++i) {
1991 gmm::copy(v.second.complex_value[i-1], v.second.complex_value[i]);
1993 gmm::copy(v.second.real_value[i-1], v.second.real_value[i]);
1994 v.second.v_num_data[i] = act_counter();
1998 bool model::is_var_newer_than_brick(
const std::string &varname,
2000 const brick_description &brick = bricks[ib];
2001 var_description &vd = variables[varname];
2002 if (niter ==
size_type(-1)) niter = vd.default_iter;
2003 return (vd.v_num > brick.v_num || vd.v_num_data[niter] > brick.v_num);
2006 bool model::is_var_mf_newer_than_brick(
const std::string &varname,
2008 const brick_description &brick = bricks[ib];
2009 var_description &vd = variables[varname];
2010 return (vd.v_num > brick.v_num);
2013 bool model::is_mim_newer_than_brick(
const mesh_im &im,
2015 const brick_description &brick = bricks[ib];
2016 return (im.version_number() > brick.v_num);
2019 void model::define_variable_group(
const std::string &group_name,
2020 const std::vector<std::string> &nl) {
2022 "variables cannot be the same as a variable name");
2024 std::set<const mesh *> ms;
2025 bool is_data_ =
false;
2026 for (
size_type i = 0; i < nl.size(); ++i) {
2031 "It is not possible to mix variables and data in a group");
2034 "All variables in a group have to exist in the model");
2036 GMM_ASSERT1(mf,
"Variables in a group should be fem variables");
2037 GMM_ASSERT1(ms.find(&(mf->linked_mesh())) == ms.end(),
2038 "Two variables in a group cannot share the same mesh");
2039 ms.insert(&(mf->linked_mesh()));
2041 variable_groups[group_name] = nl;
2044 void model::add_assembly_assignments(
const std::string &varname,
2047 GMM_ASSERT1(order < 3 || order ==
size_type(-1),
"Bad order value");
2048 const im_data *imd = pim_data_of_variable(varname);
2049 GMM_ASSERT1(imd != 0,
"Only applicable to im_data");
2050 assignement_desc as;
2051 as.varname = varname; as.expr = expr; as.region = rg; as.order = order;
2053 assignments.push_back(as);
2056 void model::add_temporaries(
const varnamelist &vl,
2057 gmm::uint64_type id_num)
const {
2058 for (
size_type i = 0; i < vl.size(); ++i) {
2059 var_description &vd = variables[vl[i]];
2060 if (vd.n_iter > 1) {
2061 vd.add_temporary(id_num);
2066 bool model::temporary_uptodate(
const std::string &varname,
2067 gmm::uint64_type id_num,
2069 var_description &vd = variables[varname];
2071 for (; ind < vd.n_iter + vd.n_temp_iter ; ++ind) {
2072 if (vd.v_num_var_iter[ind] == id_num)
break;
2074 if (ind < vd.n_iter + vd.n_temp_iter) {
2075 if (vd.v_num_iter[ind] <= vd.v_num_data[vd.default_iter]) {
2076 vd.v_num_iter[ind] = act_counter();
2085 void model::set_default_iter_of_variable(
const std::string &varname,
2088 var_description &vd = variables[varname];
2089 GMM_ASSERT1(ind < vd.n_iter + vd.n_temp_iter,
2090 "Inexistent iteration " << ind);
2091 vd.default_iter = ind;
2095 void model::reset_default_iter_of_variables(
const varnamelist &vl)
const {
2096 for (
size_type i = 0; i < vl.size(); ++i)
2097 variables[vl[i]].default_iter = 0;
2100 const model_real_sparse_matrix &
2102 GMM_ASSERT1(bricks[ib].tlist[iterm].is_matrix_term,
2103 "Not a matrix term !");
2104 GMM_ASSERT1(bricks[ib].pbr->is_linear(),
"Nonlinear term !");
2105 return bricks[ib].rmatlist[iterm];
2108 const model_complex_sparse_matrix &
2110 GMM_ASSERT1(bricks[ib].tlist[iterm].is_matrix_term,
2111 "Not a matrix term !");
2112 GMM_ASSERT1(bricks[ib].pbr->is_linear(),
"Nonlinear term !");
2113 return bricks[ib].cmatlist[iterm];
2117 void model::update_brick(
size_type ib, build_version version)
const {
2118 const brick_description &brick = bricks[ib];
2119 bool cplx =
is_complex() && brick.pbr->is_complex();
2120 bool tobecomputed = brick.terms_to_be_computed
2121 || brick.pbr->is_to_be_computed_each_time()
2122 || !(brick.pbr->is_linear());
2125 if (!tobecomputed ) {
2126 for (
size_type i = 0; i < brick.vlist.size() && !tobecomputed; ++i) {
2127 var_description &vd = variables[brick.vlist[i]];
2128 if (vd.v_num > brick.v_num) tobecomputed =
true;
2133 for (
size_type i = 0; i < brick.dlist.size() && !tobecomputed; ++i) {
2134 var_description &vd = variables[brick.dlist[i]];
2135 if (vd.v_num > brick.v_num || vd.v_num_data[vd.default_iter] > brick.v_num) {
2136 tobecomputed =
true;
2137 version = build_version(version | BUILD_ON_DATA_CHANGE);
2142 if (!tobecomputed ) {
2143 for (
size_type i = 0; i < brick.mims.size() && !tobecomputed; ++i) {
2144 if (brick.mims[i]->version_number() > brick.v_num) tobecomputed =
true;
2149 brick.external_load = scalar_type(0);
2151 if (!(brick.pdispatch))
2152 { brick_call(ib, version, 0); }
2155 brick.pdispatch->asm_complex_tangent_terms
2156 (*
this, ib, brick.cmatlist, brick.cveclist, brick.cveclist_sym,
2159 brick.pdispatch->asm_real_tangent_terms
2160 (*
this, ib, brick.rmatlist, brick.rveclist, brick.rveclist_sym,
2163 brick.v_num = act_counter();
2166 if (brick.pbr->is_linear()) brick.terms_to_be_computed =
false;
2173 const brick_description &brick = bricks[ib];
2174 if (brick.pbr->is_linear()) {
2176 bool cplx =
is_complex() && brick.pbr->is_complex();
2178 for (
size_type j = 0; j < brick.tlist.size(); ++j) {
2179 const term_description &term = brick.tlist[j];
2180 bool isg = term.is_global;
2183 size_type n_iter_1 = n_iter, n_iter_2 = n_iter;
2185 n_iter_1 = variables[term.var1].default_iter;
2186 if (term.is_matrix_term)
2187 n_iter_2 = variables[term.var2].default_iter;
2192 if (term.is_matrix_term) {
2195 model_complex_plain_vector V(nbgdof);
2196 for (VAR_SET::iterator it = variables.begin();
2197 it != variables.end(); ++it)
2198 if (it->second.is_variable) {
2200 ? it->second.default_iter : n_iter;
2201 gmm::copy(it->second.complex_value[n_iter_i],
2202 gmm::sub_vector(V, it->second.I));
2206 gmm::scaled(V, complex_type(-1)),
2207 brick.cveclist[ind_data][j]);
2211 gmm::scaled(variables[term.var2].complex_value[n_iter_2],
2213 brick.cveclist[ind_data][j]);
2217 model_real_plain_vector V(nbgdof);
2218 for (VAR_SET::iterator it = variables.begin();
2219 it != variables.end(); ++it)
2220 if (it->second.is_variable) {
2222 ? it->second.default_iter : n_iter;
2223 gmm::copy(it->second.real_value[n_iter_i],
2224 gmm::sub_vector(V, it->second.I));
2227 (brick.rmatlist[j], gmm::scaled(V, scalar_type(-1)),
2228 brick.rveclist[ind_data][j]);
2232 gmm::scaled(variables[term.var2].real_value[n_iter_2],
2233 scalar_type(-1)), brick.rveclist[ind_data][j]);
2236 if (term.is_symmetric && term.var1.compare(term.var2)) {
2239 (gmm::conjugated(brick.cmatlist[j]),
2240 gmm::scaled(variables[term.var1].complex_value[n_iter_1],
2242 brick.cveclist_sym[ind_data][j]);
2245 (gmm::transposed(brick.rmatlist[j]),
2246 gmm::scaled(variables[term.var1].real_value[n_iter_1],
2248 brick.rveclist_sym[ind_data][j]);
2255 void model::update_affine_dependent_variables() {
2256 for (VAR_SET::iterator it = variables.begin(); it != variables.end(); ++it)
2257 if (it->second.is_affine_dependent) {
2258 VAR_SET::iterator it2 = variables.find(it->second.org_name);
2259 if (it->second.size() != it2->second.size())
2260 it->second.set_size();
2261 if (it->second.is_complex) {
2262 gmm::add(gmm::scaled(it2->second.complex_value[0],
2263 complex_type(it->second.alpha)),
2264 it->second.affine_complex_value,
2265 it->second.complex_value[0]);
2267 gmm::add(gmm::scaled(it2->second.real_value[0], it->second.alpha),
2268 it->second.affine_real_value, it->second.real_value[0]);
2270 it->second.v_num = std::max(it->second.v_num, it2->second.v_num);
2271 for (
size_type i = 0; i < it->second.n_iter; ++i)
2273 it->second.v_num_data[i] = std::max(it->second.v_num_data[i],
2274 it2->second.v_num_data[i]);
2285 GMM_ASSERT1(mf,
"Works only with fem variables.");
2289 for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) {
2290 brick_description &brick = bricks[ib];
2292 bool detected =
false;
2293 for (
size_type i = 0; i < brick.vlist.size(); ++i)
2294 if (brick.vlist[i].compare(varname) == 0)
2295 { detected =
true;
break; }
2297 if (detected && brick.mims.size()) {
2299 for (
auto &pmim : brick.mims)
2302 pmim->linked_mesh()));
2303 GMM_ASSERT1(ifo >= 0,
2304 "The given region is only partially covered by "
2305 "region of brick \"" << brick.pbr->brick_name()
2306 <<
"\". Please subdivise the region");
2308 std::string expr = brick.pbr->declare_volume_assembly_string
2309 (*
this, ib, brick.vlist, brick.dlist);
2311 ga_workspace workspace(*
this);
2312 size_type order = workspace.add_expression
2313 (expr, dummy_mim, region);
2314 GMM_ASSERT1(order <= 1,
"Wrong order for a Neumann term");
2315 expr = workspace.extract_Neumann_term(varname);
2318 result +=
" + " + workspace.extract_Neumann_term(varname);
2320 result = workspace.extract_Neumann_term(varname);
2332 GMM_ASSERT1(version != BUILD_ON_DATA_CHANGE,
2333 "Invalid assembly version BUILD_ON_DATA_CHANGE");
2334 GMM_ASSERT1(version != BUILD_WITH_LIN,
2335 "Invalid assembly version BUILD_WITH_LIN");
2336 GMM_ASSERT1(version != BUILD_WITH_INTERNAL,
2337 "Invalid assembly version BUILD_WITH_INTERNAL");
2339 #if GETFEM_PARA_LEVEL > 0
2340 double t_ref = MPI_Wtime();
2342 MPI_Comm_rank(MPI_COMM_WORLD, &rk);
2343 MPI_Comm_size(MPI_COMM_WORLD, &nbp);
2346 context_check();
if (act_size_to_be_done) actualize_sizes();
2348 if (version & BUILD_MATRIX) gmm::clear(cTM);
2349 if (version & BUILD_RHS) gmm::clear(crhs);
2352 if (version & BUILD_MATRIX) gmm::clear(rTM);
2353 if (version & BUILD_RHS) gmm::clear(rrhs);
2355 clear_dof_constraints();
2356 generic_expressions.clear();
2357 update_affine_dependent_variables();
2359 if (version & BUILD_RHS) approx_external_load_ = scalar_type(0);
2361 for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) {
2363 brick_description &brick = bricks[ib];
2366 bool auto_disabled_brick =
true;
2367 for (
size_type j = 0; j < brick.vlist.size(); ++j) {
2369 auto_disabled_brick =
false;
2371 if (auto_disabled_brick)
continue;
2373 update_brick(ib, version);
2375 bool cplx =
is_complex() && brick.pbr->is_complex();
2377 scalar_type coeff0 = scalar_type(1);
2378 if (brick.pdispatch) coeff0 = brick.matrix_coeff;
2382 for (
size_type j = 0; j < brick.tlist.size(); ++j) {
2383 term_description &term = brick.tlist[j];
2384 bool isg = term.is_global, isprevious =
false;
2386 scalar_type alpha = coeff0, alpha1 = coeff0, alpha2 = coeff0;
2387 gmm::sub_interval I1(0,nbgdof), I2(0,nbgdof);
2388 var_description *var1=
nullptr, *var2=
nullptr;
2390 VAR_SET::iterator it1 = variables.find(term.var1);
2391 GMM_ASSERT1(it1 != variables.end(),
"Internal error");
2392 var1 = &(it1->second);
2393 GMM_ASSERT1(var1->is_variable,
"Assembly of data not allowed");
2395 if (term.is_matrix_term) {
2396 VAR_SET::iterator it2 = variables.find(term.var2);
2397 GMM_ASSERT1(it2 != variables.end(),
"Internal error");
2398 var2 = &(it2->second);
2400 if (!(var2->is_variable)) {
2401 std::string vorgname = sup_previous_and_dot_to_varname(term.var2);
2402 VAR_SET::iterator it3 = variables.find(vorgname);
2403 GMM_ASSERT1(it3->second.is_variable,
2404 "Assembly of data not allowed");
2408 alpha *= var1->alpha * var2->alpha;
2409 alpha1 *= var1->alpha;
2410 alpha2 *= var2->alpha;
2415 if (term.is_matrix_term && (version & BUILD_MATRIX) && !isprevious
2416 && (isg || (var1->is_enabled() && var2->is_enabled()))) {
2417 gmm::add(gmm::scaled(brick.cmatlist[j], alpha),
2418 gmm::sub_matrix(cTM, I1, I2));
2419 if (term.is_symmetric && I1.first() != I2.first())
2420 gmm::add(gmm::scaled(gmm::transposed(brick.cmatlist[j]), alpha),
2421 gmm::sub_matrix(cTM, I2, I1));
2423 if (version & BUILD_RHS) {
2425 if (isg || var1->is_enabled()) {
2426 if (brick.pdispatch)
2427 for (
size_type k = 0; k < brick.nbrhs; ++k)
2428 gmm::add(gmm::scaled(brick.cveclist[k][j],
2430 gmm::sub_vector(crhs, I1));
2432 gmm::add(gmm::scaled(brick.cveclist[0][j],
2433 complex_type(alpha1)),
2434 gmm::sub_vector(crhs, I1));
2436 if (term.is_matrix_term && brick.pbr->is_linear() &&
is_linear()) {
2437 if (var2->is_affine_dependent && var1->is_enabled())
2438 gmm::mult_add(brick.cmatlist[j],
2439 gmm::scaled(var2->affine_complex_value,
2440 complex_type(-alpha1)),
2441 gmm::sub_vector(crhs, I1));
2442 if (term.is_symmetric && I1.first() != I2.first()
2443 && var1->is_affine_dependent && var2->is_enabled())
2444 gmm::mult_add(gmm::conjugated(brick.cmatlist[j]),
2445 gmm::scaled(var1->affine_complex_value,
2446 complex_type(-alpha2)),
2447 gmm::sub_vector(crhs, I2));
2449 if (term.is_matrix_term && brick.pbr->is_linear()
2450 && (!
is_linear() || (version & BUILD_WITH_LIN))
2451 && var1->is_enabled())
2452 gmm::mult_add(brick.cmatlist[j],
2453 gmm::scaled(var2->complex_value[0],
2454 complex_type(-alpha1)),
2455 gmm::sub_vector(crhs, I1));
2456 if (term.is_symmetric && I1.first() != I2.first()
2457 && var2->is_enabled()) {
2458 if (brick.pdispatch)
2459 for (
size_type k = 0; k < brick.nbrhs; ++k)
2460 gmm::add(gmm::scaled(brick.cveclist_sym[k][j],
2462 gmm::sub_vector(crhs, I2));
2464 gmm::add(gmm::scaled(brick.cveclist_sym[0][j],
2465 complex_type(alpha2)),
2466 gmm::sub_vector(crhs, I2));
2467 if (brick.pbr->is_linear()
2468 && (!
is_linear() || (version & BUILD_WITH_LIN)))
2469 gmm::mult_add(gmm::conjugated(brick.cmatlist[j]),
2470 gmm::scaled(var1->complex_value[0],
2471 complex_type(-alpha2)),
2472 gmm::sub_vector(crhs, I2));
2476 if (term.is_matrix_term && (version & BUILD_MATRIX) && !isprevious
2477 && (isg || (var1->is_enabled() && var2->is_enabled()))) {
2478 gmm::add(gmm::scaled(brick.rmatlist[j], alpha),
2479 gmm::sub_matrix(cTM, I1, I2));
2480 if (term.is_symmetric && I1.first() != I2.first())
2481 gmm::add(gmm::scaled(gmm::transposed(brick.rmatlist[j]), alpha),
2482 gmm::sub_matrix(cTM, I2, I1));
2484 if (version & BUILD_RHS) {
2486 if (isg || var1->is_enabled()) {
2487 if (brick.pdispatch)
2488 for (
size_type k = 0; k < brick.nbrhs; ++k)
2489 gmm::add(gmm::scaled(brick.rveclist[k][j],
2491 gmm::sub_vector(crhs, I1));
2493 gmm::add(gmm::scaled(brick.rveclist[0][j], alpha1),
2494 gmm::sub_vector(crhs, I1));
2496 if (term.is_matrix_term && brick.pbr->is_linear() &&
is_linear()) {
2497 if (var2->is_affine_dependent && var1->is_enabled())
2498 gmm::mult_add(brick.rmatlist[j],
2499 gmm::scaled(var2->affine_complex_value,
2500 complex_type(-alpha1)),
2501 gmm::sub_vector(crhs, I1));
2502 if (term.is_symmetric && I1.first() != I2.first()
2503 && var1->is_affine_dependent && var2->is_enabled())
2504 gmm::mult_add(gmm::transposed(brick.rmatlist[j]),
2505 gmm::scaled(var1->affine_complex_value,
2506 complex_type(-alpha2)),
2507 gmm::sub_vector(crhs, I2));
2509 if (term.is_matrix_term && brick.pbr->is_linear()
2510 && (!
is_linear() || (version & BUILD_WITH_LIN))
2511 && var1->is_enabled())
2512 gmm::mult_add(brick.rmatlist[j],
2513 gmm::scaled(var2->complex_value[0],
2514 complex_type(-alpha1)),
2515 gmm::sub_vector(crhs, I1));
2516 if (term.is_symmetric && I1.first() != I2.first()
2517 && var2->is_enabled()) {
2518 if (brick.pdispatch)
2519 for (
size_type k = 0; k < brick.nbrhs; ++k)
2520 gmm::add(gmm::scaled(brick.rveclist_sym[k][j],
2522 gmm::sub_vector(crhs, I2));
2524 gmm::add(gmm::scaled(brick.rveclist_sym[0][j], alpha2),
2525 gmm::sub_vector(crhs, I2));
2527 if (brick.pbr->is_linear()
2528 && (!
is_linear() || (version & BUILD_WITH_LIN)))
2529 gmm::mult_add(gmm::transposed(brick.rmatlist[j]),
2530 gmm::scaled(var1->complex_value[0],
2531 complex_type(-alpha2)),
2532 gmm::sub_vector(crhs, I2));
2536 if (term.is_matrix_term && (version & BUILD_MATRIX) && !isprevious
2537 && (isg || (var1->is_enabled() && var2->is_enabled()))) {
2538 gmm::add(gmm::scaled(brick.rmatlist[j], alpha),
2539 gmm::sub_matrix(rTM, I1, I2));
2540 if (term.is_symmetric && I1.first() != I2.first())
2541 gmm::add(gmm::scaled(gmm::transposed(brick.rmatlist[j]), alpha),
2542 gmm::sub_matrix(rTM, I2, I1));
2544 if (version & BUILD_RHS) {
2546 auto vec_out1 = gmm::sub_vector(rrhs, I1);
2547 if (isg || var1->is_enabled()) {
2548 if (brick.pdispatch)
2549 for (
size_type k = 0; k < brick.nbrhs; ++k)
2550 gmm::add(gmm::scaled(brick.rveclist[k][j],
2554 gmm::add(gmm::scaled(brick.rveclist[0][j], alpha1),
2557 if (var1->is_enabled()
2558 && term.is_matrix_term && brick.pbr->is_linear()) {
2559 bool affine_contrib(
is_linear() && var2->is_affine_dependent);
2560 bool linear_contrib(!
is_linear() || (version & BUILD_WITH_LIN));
2561 const auto &matj = brick.rmatlist[j];
2562 const auto vec_affine2 = gmm::scaled(var2->affine_real_value,
2564 const auto vec_linear2 = gmm::scaled(var2->real_value[0],
2567 model_real_plain_vector vec_tmp1(I1.size(), 0.);
2569 gmm::mult(matj, vec_affine2, vec_tmp1);
2571 gmm::mult_add(matj, vec_linear2, vec_tmp1);
2572 MPI_SUM_VECTOR(vec_tmp1);
2573 gmm::add(vec_tmp1, vec_out1);
2576 gmm::mult_add(matj, vec_affine2, vec_out1);
2578 gmm::mult_add(matj, vec_linear2, vec_out1);
2583 if (term.is_symmetric && I1.first() != I2.first() &&
2584 var2->is_enabled()) {
2585 auto vec_out2 = gmm::sub_vector(rrhs, I2);
2586 if (brick.pdispatch)
2587 for (
size_type k = 0; k < brick.nbrhs; ++k)
2588 gmm::add(gmm::scaled(brick.rveclist_sym[k][j],
2592 gmm::add(gmm::scaled(brick.rveclist_sym[0][j], alpha2),
2594 if (term.is_matrix_term && brick.pbr->is_linear()) {
2595 bool affine_contrib(
is_linear() && var1->is_affine_dependent);
2596 bool linear_contrib(!
is_linear() || (version & BUILD_WITH_LIN));
2597 const auto matj_trans = gmm::transposed(brick.rmatlist[j]);
2598 const auto vec_affine1 = gmm::scaled(var1->affine_real_value,
2600 const auto vec_linear1 = gmm::scaled(var1->real_value[0],
2603 model_real_plain_vector vec_tmp2(I2.size(),0.);
2605 gmm::mult(matj_trans, vec_affine1, vec_tmp2);
2607 gmm::mult_add(matj_trans, vec_linear1, vec_tmp2);
2608 MPI_SUM_VECTOR(vec_tmp2);
2609 gmm::add(vec_tmp2, vec_out2);
2612 gmm::mult_add(matj_trans, vec_affine1, vec_out2);
2614 gmm::mult_add(matj_trans, vec_linear1, vec_out2);
2622 if (brick.pbr->is_linear())
2623 brick.terms_to_be_computed =
false;
2635 if (version & BUILD_RHS) approx_external_load_ += brick.external_load;
2638 if (version & BUILD_RHS && version & BUILD_WITH_INTERNAL) {
2641 gmm::sub_interval IP(0,gmm::vect_size(rrhs));
2642 gmm::fill(full_rrhs, 0.);
2643 gmm::copy(rrhs, gmm::sub_vector(full_rrhs, IP));
2647 if (generic_expressions.size()) {
2650 if (version & BUILD_RHS)
2651 GMM_TRACE2(
"Global generic assembly RHS");
2652 if (version & BUILD_MATRIX)
2653 GMM_TRACE2(
"Global generic assembly tangent term");
2656 auto add_assignments_and_expressions_to_workspace =
2657 [&](ga_workspace &workspace)
2659 for (
const auto &ad : assignments)
2660 workspace.add_assignment_expression
2661 (ad.varname, ad.expr, ad.region, ad.order, ad.before);
2662 for (
const auto &ge : generic_expressions)
2663 workspace.add_expression(ge.expr, ge.mim, ge.region,
2664 2, ge.secondary_domain);
2667 const bool with_internal = version & BUILD_WITH_INTERNAL
2669 model_real_sparse_matrix intern_mat;
2670 model_real_plain_vector res0,
2673 size_type full_size = gmm::vect_size(full_rrhs),
2674 primary_size = gmm::vect_size(rrhs);
2676 if ((version & BUILD_RHS) || (version & BUILD_MATRIX && with_internal))
2677 gmm::resize(res0, with_internal ? full_size : primary_size);
2678 if (version & BUILD_MATRIX && with_internal)
2679 gmm::resize(res1, full_size);
2681 if (version & BUILD_MATRIX) {
2682 if (with_internal) {
2683 gmm::resize(intern_mat, full_size, primary_size);
2684 gmm::resize(res1, full_size);
2690 if (version & BUILD_RHS) {
2693 ga_workspace workspace(*
this);
2694 add_assignments_and_expressions_to_workspace(workspace);
2695 workspace.set_assembled_vector(res0_distro);
2696 workspace.assembly(1, with_internal);
2697 if (with_internal) {
2698 gmm::copy(res0_distro.get(), res1_distro.get());
2699 gmm::add(gmm::scaled(full_rrhs, scalar_type(-1)),
2701 workspace.set_assembled_vector(res1_distro);
2702 workspace.set_internal_coupling_matrix(intern_mat_distro);
2704 workspace.set_assembled_matrix(tangent_matrix_distro);
2705 workspace.assembly(2, with_internal);
2710 ga_workspace workspace(*
this);
2711 add_assignments_and_expressions_to_workspace(workspace);
2712 if (with_internal) {
2713 gmm::copy(gmm::scaled(full_rrhs, scalar_type(-1)),
2715 workspace.set_assembled_vector(res1_distro);
2716 workspace.set_internal_coupling_matrix(intern_mat_distro);
2718 workspace.set_assembled_matrix(tangent_matrix_distro);
2719 workspace.assembly(2, with_internal);
2723 else if (version & BUILD_RHS) {
2726 ga_workspace workspace(*
this);
2727 add_assignments_and_expressions_to_workspace(workspace);
2728 workspace.set_assembled_vector(res0_distro);
2729 workspace.assembly(1, with_internal);
2733 if (version & BUILD_RHS) {
2734 gmm::scale(res0, scalar_type(-1));
2735 if (with_internal) {
2736 gmm::sub_interval IP(0,gmm::vect_size(rrhs));
2737 gmm::add(gmm::sub_vector(res0, IP), rrhs);
2738 gmm::add(res0, full_rrhs);
2740 gmm::add(res0, rrhs);
2743 if (version & BUILD_MATRIX && with_internal) {
2744 gmm::scale(res1, scalar_type(-1));
2745 gmm::sub_interval IP(0, primary_size),
2746 II(primary_size, full_size-primary_size);
2747 gmm::copy(gmm::sub_matrix(intern_mat, II, IP), internal_rTM);
2748 gmm::add(gmm::sub_vector(res1, IP), rrhs);
2749 gmm::copy(gmm::sub_vector(res1, II), internal_sol);
2754 if ((version & BUILD_RHS) || (version & BUILD_MATRIX)) {
2756 std::vector<size_type> dof_indices;
2757 std::vector<complex_type> dof_pr_values;
2758 std::vector<complex_type> dof_go_values;
2759 for (
const auto &keyval : complex_dof_constraints) {
2760 const gmm::sub_interval &I = interval_of_variable(keyval.first);
2762 for (
const auto &val : keyval.second) {
2763 dof_indices.push_back(val.first + I.first());
2764 dof_go_values.push_back(val.second);
2765 dof_pr_values.push_back(V[val.first]);
2769 if (dof_indices.size()) {
2770 gmm::sub_index SI(dof_indices);
2771 gmm::sub_interval II(0,
nb_dof());
2773 if (version & BUILD_RHS) {
2774 if (MPI_IS_MASTER())
2775 approx_external_load_ += gmm::vect_norm1(dof_go_values);
2777 if (is_symmetric_) {
2778 scalar_type valnorm = gmm::vect_norm2(dof_go_values);
2779 if (valnorm > scalar_type(0)) {
2780 GMM_ASSERT1(version & BUILD_MATRIX,
"Rhs only for a "
2781 "symmetric linear problem with dof "
2782 "constraint not allowed");
2783 model_complex_plain_vector vv(gmm::vect_size(crhs));
2784 gmm::mult(gmm::sub_matrix(cTM, II, SI), dof_go_values, vv);
2786 gmm::add(gmm::scaled(vv, scalar_type(-1)), crhs);
2789 gmm::copy(dof_go_values, gmm::sub_vector(crhs, SI));
2791 gmm::add(dof_go_values,
2792 gmm::scaled(dof_pr_values, complex_type(-1)),
2793 gmm::sub_vector(crhs, SI));
2796 if (version & BUILD_MATRIX) {
2797 gmm::clear(gmm::sub_matrix(cTM, SI, II));
2798 if (is_symmetric_) gmm::clear(gmm::sub_matrix(cTM, II, SI));
2800 if (MPI_IS_MASTER()) {
2801 for (
size_type i = 0; i < dof_indices.size(); ++i)
2802 cTM(dof_indices[i], dof_indices[i]) = complex_type(1);
2807 std::vector<size_type> dof_indices;
2808 std::vector<scalar_type> dof_pr_values;
2809 std::vector<scalar_type> dof_go_values;
2810 for (
const auto &keyval : real_dof_constraints) {
2811 const gmm::sub_interval &I = interval_of_variable(keyval.first);
2812 const model_real_plain_vector &V =
real_variable(keyval.first);
2813 for (
const auto &val : keyval.second) {
2814 dof_indices.push_back(val.first + I.first());
2815 dof_go_values.push_back(val.second);
2816 dof_pr_values.push_back(V[val.first]);
2820 #if GETFEM_PARA_LEVEL > 1
2821 GMM_ASSERT1(MPI_IS_MASTER() || (dof_indices.size() == 0),
2822 "Sorry, for the moment, the dof constraints have to be "
2823 "added on the master process only");
2824 size_type dof_indices_size = dof_indices.size();
2825 MPI_BCAST0_SCALAR(dof_indices_size);
2826 dof_indices.resize(dof_indices_size);
2827 MPI_BCAST0_VECTOR(dof_indices);
2828 dof_pr_values.resize(dof_indices_size);
2829 MPI_BCAST0_VECTOR(dof_pr_values);
2830 dof_go_values.resize(dof_indices_size);
2831 MPI_BCAST0_VECTOR(dof_go_values);
2834 if (dof_indices.size()) {
2835 gmm::sub_index SI(dof_indices);
2836 gmm::sub_interval II(0,
nb_dof());
2838 if (version & BUILD_RHS) {
2839 if (MPI_IS_MASTER())
2840 approx_external_load_ += gmm::vect_norm1(dof_go_values);
2842 if (is_symmetric_) {
2843 scalar_type valnorm = gmm::vect_norm2(dof_go_values);
2844 if (valnorm > scalar_type(0)) {
2845 GMM_ASSERT1(version & BUILD_MATRIX,
"Rhs only for a "
2846 "symmetric linear problem with dof "
2847 "constraint not allowed");
2848 model_real_plain_vector vv(gmm::vect_size(rrhs));
2849 gmm::mult(gmm::sub_matrix(rTM, II, SI), dof_go_values, vv);
2851 gmm::add(gmm::scaled(vv, scalar_type(-1)), rrhs);
2854 gmm::copy(dof_go_values, gmm::sub_vector(rrhs, SI));
2856 gmm::add(dof_go_values,
2857 gmm::scaled(dof_pr_values, scalar_type(-1)),
2858 gmm::sub_vector(rrhs, SI));
2861 if (version & BUILD_MATRIX) {
2862 gmm::clear(gmm::sub_matrix(rTM, SI, II));
2863 if (is_symmetric_) gmm::clear(gmm::sub_matrix(rTM, II, SI));
2865 if (MPI_IS_MASTER()) {
2866 for (
size_type i = 0; i < dof_indices.size(); ++i)
2867 rTM(dof_indices[i], dof_indices[i]) = scalar_type(1);
2874 if (version & BUILD_RHS) {
2877 MPI_BCAST0_SCALAR(approx_external_load_);
2880 #if GETFEM_PARA_LEVEL > 0
2881 if (MPI_IS_MASTER()) cout <<
"Assembly time " << MPI_Wtime()-t_ref << endl;
2890 return it->second.associated_mf();
2896 return it->second.passociated_mf();
2900 model::qdims_of_variable(
const std::string &name)
const {
2902 const mesh_fem *mf = it->second.passociated_mf();
2903 const im_data *imd = it->second.imd;
2906 bgeot::multi_index mi = mf->get_qdims();
2907 if (n > 1 || it->second.qdims.size() > 1) {
2909 if (mi.back() == 1) { mi.back() *= it->second.qdims[0]; ++i; }
2910 for (; i < it->second.qdims.size(); ++i)
2911 mi.push_back(it->second.qdims[i]);
2915 bgeot::multi_index mi = imd->tensor_size();
2916 if (n > 1 || it->second.qdims.size() > 1) {
2918 if (mi.back() == 1) { mi.back() *= it->second.qdims[0]; ++i; }
2919 for (; i < it->second.qdims.size(); ++i)
2920 mi.push_back(it->second.qdims[i]);
2924 return it->second.qdims;
2927 size_type model::qdim_of_variable(
const std::string &name)
const {
2929 const mesh_fem *mf = it->second.passociated_mf();
2930 const im_data *imd = it->second.imd;
2933 return mf->get_qdim() * n;
2935 return imd->tensor_size().total_size() * n;
2941 const gmm::sub_interval &
2942 model::interval_of_variable(
const std::string &name)
const {
2944 if (act_size_to_be_done) actualize_sizes();
2945 VAR_SET::const_iterator it = find_variable(name);
2946 return it->second.I;
2949 const model_real_plain_vector &
2955 const model_real_plain_vector &
2957 GMM_ASSERT1(!complex_version,
"This model is a complex one");
2958 GMM_ASSERT1(!
is_old(name),
"Please don't use Old_ prefix in combination "
2959 "with variable version");
2961 auto it = variables.find(name);
2962 GMM_ASSERT1(it != variables.end(),
"Undefined variable " << name);
2963 if (act_size_to_be_done && it->second.mf) {
2964 if (it->second.filter != VDESCRFILTER_NO)
2967 it->second.set_size();
2969 if (niter ==
size_type(-1)) niter = it->second.default_iter;
2970 GMM_ASSERT1(it->second.n_iter + it->second.n_temp_iter > niter,
2971 "Invalid iteration number " << niter <<
" for " << name);
2972 return it->second.real_value[niter];
2975 const model_complex_plain_vector &
2981 const model_complex_plain_vector &
2983 GMM_ASSERT1(complex_version,
"This model is a real one");
2984 GMM_ASSERT1(!
is_old(name),
"Please don't use Old_ prefix in combination with"
2985 " variable version");
2987 auto it = variables.find(name);
2988 GMM_ASSERT1(it!=variables.end(),
"Undefined variable " << name);
2989 if (act_size_to_be_done && it->second.mf) {
2990 if (it->second.filter != VDESCRFILTER_NO)
2993 it->second.set_size();
2995 if (niter ==
size_type(-1)) niter = it->second.default_iter;
2996 GMM_ASSERT1(it->second.n_iter + it->second.n_temp_iter > niter,
2997 "Invalid iteration number "
2998 << niter <<
" for " << name);
2999 return it->second.complex_value[niter];
3002 model_real_plain_vector &
3009 model_real_plain_vector &
3011 GMM_ASSERT1(!complex_version,
"This model is a complex one");
3012 GMM_ASSERT1(!
is_old(name),
"Please don't use Old_ prefix in combination with"
3013 " variable version");
3015 auto it = variables.find(name);
3016 GMM_ASSERT1(it!=variables.end(),
"Undefined variable " << name);
3017 if (act_size_to_be_done && it->second.mf) {
3018 if (it->second.filter != VDESCRFILTER_NO)
3021 it->second.set_size();
3023 if (niter ==
size_type(-1)) niter = it->second.default_iter;
3024 it->second.v_num_data[niter] = act_counter();
3025 GMM_ASSERT1(it->second.n_iter + it->second.n_temp_iter > niter,
3026 "Invalid iteration number "
3027 << niter <<
" for " << name);
3028 return it->second.real_value[niter];
3031 model_complex_plain_vector &
3037 model_complex_plain_vector &
3039 GMM_ASSERT1(complex_version,
"This model is a real one");
3040 GMM_ASSERT1(!
is_old(name),
"Please don't use Old_ prefix in combination with"
3041 " variable version");
3043 auto it = variables.find(name);
3044 GMM_ASSERT1(it!=variables.end(),
"Undefined variable " << name);
3045 if (act_size_to_be_done && it->second.mf) {
3046 if (it->second.filter != VDESCRFILTER_NO)
3049 it->second.set_size();
3051 if (niter ==
size_type(-1)) niter = it->second.default_iter;
3052 it->second.v_num_data[niter] = act_counter();
3053 GMM_ASSERT1(it->second.n_iter + it->second.n_temp_iter > niter,
3054 "Invalid iteration number "
3055 << niter <<
" for " << name);
3056 return it->second.complex_value[niter];
3059 model_real_plain_vector &
3060 model::set_real_constant_part(
const std::string &name)
const {
3061 GMM_ASSERT1(!complex_version,
"This model is a complex one");
3063 VAR_SET::iterator it = variables.find(name);
3064 GMM_ASSERT1(it != variables.end(),
"Undefined variable " << name);
3065 GMM_ASSERT1(it->second.is_affine_dependent,
3066 "Only for affine dependent variables");
3067 if (act_size_to_be_done && it->second.mf) {
3068 if (it->second.filter != VDESCRFILTER_NO)
3071 it->second.set_size();
3073 for (
auto &v_num : it->second.v_num_data) v_num = act_counter();
3074 return it->second.affine_real_value;
3077 model_complex_plain_vector &
3078 model::set_complex_constant_part(
const std::string &name)
const {
3079 GMM_ASSERT1(complex_version,
"This model is a real one");
3081 VAR_SET::iterator it = variables.find(name);
3082 GMM_ASSERT1(it!=variables.end(),
"Undefined variable " << name);
3083 if (act_size_to_be_done && it->second.mf) {
3084 if (it->second.filter != VDESCRFILTER_NO)
3087 it->second.set_size();
3089 for (
auto &v_num : it->second.v_num_data) v_num = act_counter();
3090 return it->second.affine_complex_value;
3095 const brick_description &brick = bricks[ind_brick];
3096 update_brick(ind_brick, model::BUILD_ALL);
3098 brick.pbr->check_stiffness_matrix_and_rhs(*
this, ind_brick, brick.tlist,
3099 brick.vlist, brick.dlist, brick.mims, brick.rmatlist,
3100 brick.rveclist[0], brick.rveclist_sym[0], brick.region);
3103 void model::clear() {
3105 active_bricks.clear();
3106 valid_bricks.clear();
3107 real_dof_constraints.clear();
3108 complex_dof_constraints.clear();
3110 rTM = model_real_sparse_matrix();
3111 cTM = model_complex_sparse_matrix();
3112 rrhs = model_real_plain_vector();
3113 crhs = model_complex_plain_vector();
3118 void virtual_brick::full_asm_real_tangent_terms_(
const model &md,
size_type ind_brick,
3119 const model::varnamelist &term_list,
3120 const model::varnamelist &var_list,
3121 const model::mimlist &mim_list,
3122 model::real_matlist &rmatlist,
3123 model::real_veclist &rveclist,
3124 model::real_veclist &rveclist_sym,
3125 size_type region, build_version version)
const
3128 rveclist, rveclist_sym, region, version);
3130 rveclist, rveclist_sym, region, version);
3132 rveclist, rveclist_sym, region, version);
3137 const model::termlist& tlist,
3138 const model::varnamelist &vl,
3139 const model::varnamelist &dl,
3140 const model::mimlist &mims,
3141 model::real_matlist &matl,
3142 model::real_veclist &rvc1,
3143 model::real_veclist &rvc2,
3145 const scalar_type TINY)
const
3147 std::cout<<
"******Verifying stiffnesses of *******"<<std::endl;
3148 std::cout<<
"*** "<<brick_name()<<std::endl;
3151 std::map<std::string,size_type> rhs_index;
3152 for(
size_type iterm=0;iterm<matl.size();iterm++)
3153 if (tlist[iterm].var1==tlist[iterm].var2) rhs_index[tlist[iterm].var1]=iterm;
3155 if (rhs_index.size()==0){
3156 GMM_WARNING0(
"*** cannot verify stiffness for this brick***");
3159 full_asm_real_tangent_terms_(md, s, vl, dl, mims, matl, rvc1, rvc2,
3160 rg, model::BUILD_MATRIX);
3161 for(
size_type iterm=0;iterm<matl.size();iterm++)
3164 std::cout<<std::endl;
3165 std::cout<<
" Stiffness["<<tlist[iterm].var1
3166 <<
","<<tlist[iterm].var2<<
"]:"<<std::endl;
3169 std::cout<<
" "<<tlist[iterm].var1<<
" has zero size. Skipping this term"<<std::endl;
3174 std::cout<<
" "<<tlist[iterm].var2<<
" has zero size. Skipping this term"<<std::endl;
3178 model_real_sparse_matrix SM(matl[iterm]);
3179 gmm::fill(rvc1[rhs_index[tlist[iterm].var1]], 0.0);
3180 full_asm_real_tangent_terms_(md, s, vl, dl, mims, matl, rvc1, rvc2,
3181 rg, model::BUILD_RHS);
3182 if (gmm::mat_euclidean_norm(matl[iterm])<1e-12){
3183 std::cout<<
" The assembled matrix is nearly zero, skipping."<<std::endl;
3186 model_real_plain_vector RHS0(rvc1[rhs_index[tlist[iterm].var1]]);
3189 model_real_sparse_matrix fdSM(matl[iterm].nrows(), matl[iterm].ncols());
3191 model_real_plain_vector& RHS1 = rvc1[rhs_index[tlist[iterm].var1]];
3193 scalar_type relative_tiny = gmm::vect_norminf(RHS1)*TINY;
3194 if (relative_tiny < TINY) relative_tiny = TINY;
3196 for (
size_type j = 0; j < matl[iterm].ncols(); j++){
3197 U[j] += relative_tiny;
3198 gmm::fill(RHS1, 0.0);
3199 full_asm_real_tangent_terms_(md, s, vl, dl, mims, matl, rvc1, rvc2,
3200 rg, model::BUILD_RHS);
3201 for (
size_type i = 0; i<matl[iterm].nrows(); i++)
3202 fdSM(i, j) = (RHS0[i] - RHS1[i]) / relative_tiny;
3203 U[j] -= relative_tiny;
3206 model_real_sparse_matrix diffSM(matl[iterm].nrows(),matl[iterm].ncols());
3207 gmm::add(SM,gmm::scaled(fdSM,-1.0),diffSM);
3208 scalar_type norm_error_euc
3209 = gmm::mat_euclidean_norm(diffSM)/gmm::mat_euclidean_norm(SM)*100;
3210 scalar_type norm_error_1
3211 = gmm::mat_norm1(diffSM)/gmm::mat_norm1(SM)*100;
3212 scalar_type norm_error_max
3213 = gmm::mat_maxnorm(diffSM)/gmm::mat_maxnorm(SM)*100;
3216 scalar_type nsym_norm_error_euc=0.0;
3217 scalar_type nsym_norm_error_1=0.0;
3218 scalar_type nsym_norm_error_max=0.0;
3219 if (tlist[iterm].var1==tlist[iterm].var2){
3220 model_real_sparse_matrix diffSMtransposed(matl[iterm].nrows(),matl[iterm].ncols());
3221 gmm::add(gmm::transposed(fdSM),gmm::scaled(fdSM,-1.0),diffSMtransposed);
3223 = gmm::mat_euclidean_norm(diffSMtransposed)/gmm::mat_euclidean_norm(fdSM)*100;
3225 = gmm::mat_norm1(diffSMtransposed)/gmm::mat_norm1(fdSM)*100;
3227 = gmm::mat_maxnorm(diffSMtransposed)/gmm::mat_maxnorm(fdSM)*100;
3231 if(rvc1[0].size()<8){
3232 std::cout <<
"RHS Stiffness Matrix: \n";
3233 std::cout <<
"------------------------\n";
3234 for(
size_type i=0; i < rvc1[iterm].size(); ++i){
3236 for(
size_type j=0; j < rvc1[iterm].size(); ++j){
3237 std::cout << fdSM(i,j) <<
" ";
3241 std::cout <<
"Analytical Stiffness Matrix: \n";
3242 std::cout <<
"------------------------\n";
3243 for(
size_type i=0; i < rvc1[iterm].size(); ++i){
3245 for(
size_type j=0; j < rvc1[iterm].size(); ++j){
3246 std::cout << matl[iterm](i,j) <<
" ";
3250 std::cout <<
"Vector U: \n";
3251 std::cout <<
"------------------------\n";
3252 for(
size_type i=0; i < rvc1[iterm].size(); ++i){
3259 <<
"\n\nfinite diff test error_norm_eucl: " << norm_error_euc <<
"%\n"
3260 <<
"finite diff test error_norm1: " << norm_error_1 <<
"%\n"
3261 <<
"finite diff test error_max_norm: " << norm_error_max <<
"%\n\n\n";
3263 if (tlist[iterm].var1==tlist[iterm].var2){
3265 <<
"Nonsymmetrical test error_norm_eucl: "<< nsym_norm_error_euc<<
"%\n"
3266 <<
"Nonsymmetrical test error_norm1: " << nsym_norm_error_1 <<
"%\n"
3267 <<
"Nonsymmetrical test error_max_norm: " << nsym_norm_error_max <<
"%"
3287 struct gen_source_term_assembly_brick :
public virtual_brick {
3289 std::string expr, directvarname, directdataname;
3290 model::varnamelist vl_test1;
3291 std::string secondary_domain;
3294 const model::varnamelist &,
3295 const model::varnamelist &,
3296 const model::mimlist &mims,
3297 model::real_matlist &,
3298 model::real_veclist &vecl,
3299 model::real_veclist &,
3301 build_version)
const override {
3302 GMM_ASSERT1(vecl.size() == vl_test1.size()
3303 + ((directdataname.size() == 0) ? 0 : 1),
"Wrong number "
3304 "of terms for Generic source term assembly brick ");
3305 GMM_ASSERT1(mims.size() == 1,
"Generic source term assembly brick "
3306 "needs one and only one mesh_im");
3307 GMM_TRACE2(
"Generic source term assembly");
3309 gmm::clear(vecl[0]);
3312 const mesh_im &mim = *mims[0];
3317 ga_workspace workspace(md, ga_workspace::inherit::ALL);
3318 GMM_TRACE2(name <<
": generic source term assembly");
3319 workspace.add_expression(expr, mim, rg, 1, secondary_domain);
3320 workspace.assembly(1);
3321 const auto &V=workspace.assembled_vector();
3322 for (
size_type i = 0; i < vl_test1.size(); ++i) {
3323 const auto &I=workspace.interval_of_variable(vl_test1[i]);
3324 gmm::copy(gmm::sub_vector(V, I), vecl[i]);
3328 if (directvarname.size()) {
3333 void real_post_assembly_in_serial(
const model &md,
size_type ib,
3334 const model::varnamelist &,
3335 const model::varnamelist &,
3336 const model::mimlist &,
3337 model::real_matlist &,
3338 model::real_veclist &vecl,
3339 model::real_veclist &,
3341 build_version)
const override {
3342 scalar_type el = scalar_type(0);
3343 for (
size_type i=0; i < vecl.size(); ++i) el += gmm::vect_norm1(vecl[i]);
3344 md.add_external_load(ib, el);
3347 virtual std::string declare_volume_assembly_string
3348 (
const model &,
size_type,
const model::varnamelist &,
3349 const model::varnamelist &)
const {
3350 return std::string();
3353 gen_source_term_assembly_brick(
const std::string &expr_,
3354 std::string brickname,
3355 const model::varnamelist &vl_test1_,
3356 const std::string &directvarname_,
3357 const std::string &directdataname_,
3358 const std::string &secdom)
3359 : vl_test1(vl_test1_), secondary_domain(secdom) {
3360 if (brickname.size() == 0)
3361 brickname =
"Generic source term assembly brick";
3363 set_flags(brickname,
true ,
3367 directvarname = directvarname_; directdataname = directdataname_;
3372 static bool check_compatibility_vl_test(model &md,
3373 const model::varnamelist vl_test) {
3374 model::varnamelist org;
3375 for (
size_type i = 0; i < vl_test.size(); ++i) {
3376 if (md.is_affine_dependent_variable(vl_test[i]))
3377 org.push_back(md.org_variable(vl_test[i]));
3379 for (
size_type i = 0; i < vl_test.size(); ++i)
3380 for (
size_type j = 0; j < org.size(); ++j)
3381 if (vl_test[i].compare(org[j]) == 0)
return false;
3386 (model &md,
const mesh_im &mim,
const std::string &expr,
size_type region,
3387 const std::string &brickname, std::string directvarname,
3388 const std::string &directdataname,
bool return_if_nonlin,
3389 const std::string &secondary_domain) {
3391 ga_workspace workspace(md);
3392 size_type order = workspace.add_expression(expr, mim, region, 1,
3394 GMM_ASSERT1(order <= 1,
"Wrong order for a source term");
3395 model::varnamelist vl, vl_test1, vl_test2, dl;
3396 bool is_lin = workspace.used_variables(vl, vl_test1, vl_test2, dl, 1);
3397 if (!is_lin && return_if_nonlin)
return size_type(-1);
3398 GMM_ASSERT1(is_lin,
"Nonlinear term");
3399 GMM_ASSERT1(check_compatibility_vl_test(md, vl_test1),
3400 "This brick do not support the assembly on both an affine "
3401 "dependent variable and its original variable. "
3402 "Split the brick.");
3404 if (directdataname.size()) {
3405 vl.push_back(directvarname);
3406 dl.push_back(directdataname);
3407 }
else directvarname =
"";
3409 pbrick pbr = std::make_shared<gen_source_term_assembly_brick>
3410 (expr, brickname, vl_test1, directvarname, directdataname,
3414 for (
size_type i = 0; i < vl_test1.size(); ++i)
3415 tl.push_back(model::term_description(vl_test1[i]));
3416 if (directdataname.size())
3417 tl.push_back(model::term_description(directvarname));
3419 return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
3424 const std::string &brickname,
const std::string &directvarname,
3425 const std::string &directdataname,
bool return_if_nonlin) {
3426 return add_source_term_(md, mim, expr, region, brickname, directvarname,
3427 directdataname, return_if_nonlin,
"");
3432 const std::string &secondary_domain,
3433 const std::string &brickname,
const std::string &directvarname,
3434 const std::string &directdataname,
bool return_if_nonlin) {
3435 return add_source_term_(md, mim, expr, region, brickname, directvarname,
3436 directdataname, return_if_nonlin, secondary_domain);
3445 struct gen_linear_assembly_brick :
public virtual_brick {
3449 model::varnamelist vl_test1, vl_test2;
3450 std::string secondary_domain;
3452 virtual void asm_real_tangent_terms(
const model &md,
size_type ib,
3453 const model::varnamelist &,
3454 const model::varnamelist &dl,
3455 const model::mimlist &mims,
3456 model::real_matlist &matl,
3457 model::real_veclist &,
3458 model::real_veclist &,
3460 build_version version)
const {
3461 GMM_ASSERT1(matl.size() == vl_test1.size(),
3462 "Wrong number of terms for Generic linear assembly brick");
3463 GMM_ASSERT1(mims.size() == 1,
3464 "Generic linear assembly brick needs one and only one "
3466 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0);
3467 for (
size_type i = 0; i < dl.size(); ++i)
3468 recompute_matrix = recompute_matrix ||
3469 md.is_var_newer_than_brick(dl[i], ib);
3471 if (recompute_matrix) {
3473 ga_workspace workspace(md, ga_workspace::inherit::ALL);
3474 workspace.add_expression(expr, *(mims[0]), region, 2, secondary_domain);
3475 GMM_TRACE2(name <<
": generic matrix assembly");
3476 workspace.assembly(2);
3477 const auto &R=workspace.assembled_matrix();
3478 for (
size_type i = 0; i < vl_test1.size(); ++i) {
3479 scalar_type alpha = scalar_type(1)
3480 / ( workspace.factor_of_variable(vl_test1[i]) *
3481 workspace.factor_of_variable(vl_test2[i]));
3482 const auto &I1=workspace.interval_of_variable(vl_test1[i]),
3483 &I2=workspace.interval_of_variable(vl_test2[i]);
3484 gmm::copy(gmm::scaled(gmm::sub_matrix(R, I1, I2), alpha),
3490 virtual std::string declare_volume_assembly_string
3491 (
const model &,
size_type,
const model::varnamelist &,
3492 const model::varnamelist &)
const {
3493 return is_lower_dim ? std::string() : expr;
3496 gen_linear_assembly_brick(
const std::string &expr_,
const mesh_im &mim,
3498 bool is_coer, std::string brickname,
3499 const model::varnamelist &vl_test1_,
3500 const model::varnamelist &vl_test2_,
3501 const std::string &secdom)
3502 : vl_test1(vl_test1_), vl_test2(vl_test2_), secondary_domain(secdom) {
3503 if (brickname.size() == 0) brickname =
"Generic linear assembly brick";
3505 is_lower_dim = mim.is_lower_dimensional();
3506 set_flags(brickname,
true ,
3513 static bool check_compatibility_vl_test(model &md,
3514 const model::varnamelist vl_test1,
3515 const model::varnamelist vl_test2) {
3516 for (
size_type i = 0; i < vl_test1.size(); ++i)
3517 for (
size_type j = 0; j < vl_test1.size(); ++j) {
3518 bool is1 = md.is_affine_dependent_variable(vl_test1[i]);
3519 bool is2 = md.is_affine_dependent_variable(vl_test2[i]);
3521 const std::string &org1
3522 = is1 ? md.org_variable(vl_test1[i]) : vl_test1[i];
3523 const std::string &org2
3524 = is2 ? md.org_variable(vl_test2[i]) : vl_test2[i];
3525 bool is1_bis = md.is_affine_dependent_variable(vl_test1[j]);
3526 bool is2_bis = md.is_affine_dependent_variable(vl_test2[j]);
3527 const std::string &org1_bis = is1_bis ? md.org_variable(vl_test1[j])
3529 const std::string &org2_bis = is2_bis ? md.org_variable(vl_test2[j])
3531 if (org1.compare(org1_bis) == 0 && org2.compare(org2_bis))
3541 (model &md,
const mesh_im &mim,
const std::string &expr,
size_type region,
3542 bool is_sym,
bool is_coercive,
const std::string &brickname,
3543 bool return_if_nonlin,
const std::string &secondary_domain) {
3545 ga_workspace workspace(md, ga_workspace::inherit::ALL);
3546 size_type order = workspace.add_expression(expr, mim, region,
3547 2, secondary_domain);
3548 model::varnamelist vl, vl_test1, vl_test2, dl;
3549 bool is_lin = workspace.used_variables(vl, vl_test1, vl_test2, dl, 2);
3551 if (!is_lin && return_if_nonlin)
return size_type(-1);
3552 GMM_ASSERT1(is_lin,
"Nonlinear term");
3553 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())
3558 brick_id = add_source_term_(md, mim, const_expr, region,
3559 brickname+
" (source term)",
3560 "",
"",
false, secondary_domain);
3564 GMM_ASSERT1(check_compatibility_vl_test(md, vl_test1, vl_test2),
3565 "This brick do not support the assembly on both an affine "
3566 "dependent variable and its original variable. "
3567 "Split the brick.");
3569 if (vl_test1.size()) {
3570 pbrick pbr = std::make_shared<gen_linear_assembly_brick>
3571 (expr, mim, is_sym, is_coercive, brickname, vl_test1, vl_test2,
3574 for (
size_type i = 0; i < vl_test1.size(); ++i)
3575 tl.push_back(model::term_description(vl_test1[i], vl_test2[i],
false));
3577 return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
3584 bool is_sym,
bool is_coercive,
const std::string &brickname,
3585 bool return_if_nonlin) {
3586 return add_linear_term_(md, mim, expr, region, is_sym, is_coercive,
3587 brickname, return_if_nonlin,
"");
3592 const std::string &secondary_domain,
bool is_sym,
bool is_coercive,
3593 const std::string &brickname,
bool return_if_nonlin) {
3594 return add_linear_term_(md, mim, expr, region, is_sym, is_coercive,
3595 brickname, return_if_nonlin, secondary_domain);
3605 struct gen_nonlinear_assembly_brick :
public virtual_brick {
3609 std::string secondary_domain;
3611 virtual void real_post_assembly_in_serial(
const model &md,
size_type ,
3612 const model::varnamelist &,
3613 const model::varnamelist &,
3614 const model::mimlist &mims,
3615 model::real_matlist &,
3616 model::real_veclist &,
3617 model::real_veclist &,
3619 build_version)
const {
3620 GMM_ASSERT1(mims.size() == 1,
3621 "Generic linear assembly brick needs one and only one "
3623 md.add_generic_expression(expr, *(mims[0]), region, secondary_domain);
3626 virtual std::string declare_volume_assembly_string
3627 (
const model &,
size_type,
const model::varnamelist &,
3628 const model::varnamelist &)
const {
3633 gen_nonlinear_assembly_brick(
const std::string &expr_,
const mesh_im &mim,
3636 std::string brickname,
3637 const std::string &secdom) {
3638 if (brickname.size() == 0) brickname =
"Generic linear assembly brick";
3640 secondary_domain = secdom;
3641 is_lower_dim = mim.is_lower_dimensional();
3642 set_flags(brickname,
false ,
3650 (model &md,
const mesh_im &mim,
const std::string &expr,
size_type region,
3651 bool is_sym,
bool is_coercive,
const std::string &brickname,
3652 const std::string &secondary_domain) {
3654 ga_workspace workspace(md);
3655 size_type order = workspace.add_expression(expr, mim, region, 2,
3657 GMM_ASSERT1(order < 2,
"Order two test functions (Test2) are not allowed"
3658 " in assembly string for nonlinear terms");
3659 model::varnamelist vl, vl_test1, vl_test2, ddl, dl;
3660 workspace.used_variables(vl, vl_test1, vl_test2, ddl, order);
3661 for (
size_type i = 0; i < ddl.size(); ++i)
3662 if (md.is_true_data(ddl[i])) dl.push_back(ddl[i]);
3663 else vl.push_back(ddl[i]);
3664 if (order == 0) { is_coercive = is_sym =
true; }
3665 pbrick pbr = std::make_shared<gen_nonlinear_assembly_brick>
3666 (expr, mim, is_sym, is_coercive, brickname, secondary_domain);
3670 return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
3676 bool is_sym,
bool is_coercive,
const std::string &brickname) {
3677 return add_nonlinear_term_(md, mim, expr, region, is_sym, is_coercive,
3683 const std::string &secondary_domain,
bool is_sym,
bool is_coercive,
3684 const std::string &brickname) {
3685 return add_nonlinear_term_(md, mim, expr, region, is_sym, is_coercive,
3686 brickname, secondary_domain);
3697 struct generic_elliptic_brick :
public virtual_brick {
3699 virtual void asm_real_tangent_terms(
const model &md,
size_type ,
3700 const model::varnamelist &vl,
3701 const model::varnamelist &dl,
3702 const model::mimlist &mims,
3703 model::real_matlist &matl,
3704 model::real_veclist &,
3705 model::real_veclist &,
3707 build_version)
const {
3708 GMM_ASSERT1(matl.size() == 1,
3709 "Generic elliptic brick has one and only one term");
3710 GMM_ASSERT1(mims.size() == 1,
3711 "Generic elliptic brick need one and only one mesh_im");
3712 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
3713 "Wrong number of variables for generic elliptic brick");
3715 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
3716 const mesh &m = mf_u.linked_mesh();
3717 size_type N = m.dim(), Q = mf_u.get_qdim(), s = 1;
3718 const mesh_im &mim = *mims[0];
3719 const model_real_plain_vector *A = 0;
3720 const mesh_fem *mf_a = 0;
3721 mesh_region rg(region);
3722 m.intersect_with_mpi_region(rg);
3723 if (dl.size() > 0) {
3724 A = &(md.real_variable(dl[0]));
3725 mf_a = md.pmesh_fem_of_variable(dl[0]);
3726 s = gmm::vect_size(*A);
3727 if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
3731 GMM_TRACE2(
"Generic elliptic term assembly");
3736 (matl[0], mim, mf_u, *mf_a, *A, rg);
3739 (matl[0], mim, mf_u, *mf_a, *A, rg);
3744 (matl[0], mim, mf_u, rg);
3747 (matl[0], mim, mf_u, rg);
3748 if (A) gmm::scale(matl[0], (*A)[0]);
3750 }
else if (s == N*N) {
3754 (matl[0], mim, mf_u, *mf_a, *A, rg);
3757 (matl[0], mim, mf_u, *mf_a, *A, rg);
3761 (matl[0], mim, mf_u, *A, rg);
3764 (matl[0], mim, mf_u, *A, rg);
3766 }
else if (s == N*N*Q*Q) {
3769 (matl[0], mim, mf_u, *mf_a, *A, rg);
3772 (matl[0], mim, mf_u, *A, rg);
3774 GMM_ASSERT1(
false,
"Bad format generic elliptic brick coefficient");
3778 virtual void real_post_assembly_in_serial(
const model &md,
size_type,
3779 const model::varnamelist &,
3780 const model::varnamelist &dl,
3781 const model::mimlist &,
3782 model::real_matlist &,
3783 model::real_veclist &,
3784 model::real_veclist &,
3786 build_version)
const
3788 const model_real_plain_vector *A = 0;
3789 const mesh_fem *mf_a = 0;
3792 A = &(md.real_variable(dl[0]));
3793 mf_a = md.pmesh_fem_of_variable(dl[0]);
3797 virtual void complex_post_assembly_in_serial(
const model &md,
size_type,
3798 const model::varnamelist &,
3799 const model::varnamelist &dl,
3800 const model::mimlist &,
3801 model::complex_matlist &,
3802 model::complex_veclist &,
3803 model::complex_veclist &,
3805 build_version)
const
3807 const model_real_plain_vector *A = 0;
3808 const mesh_fem *mf_a = 0;
3811 A = &(md.real_variable(dl[0]));
3812 mf_a = md.pmesh_fem_of_variable(dl[0]);
3816 virtual scalar_type asm_complex_tangent_terms(
const model &md,
size_type,
3817 const model::varnamelist &vl,
3818 const model::varnamelist &,
3819 const model::mimlist &,
3820 model::complex_matlist &matl,
3821 model::complex_veclist &,
3822 model::complex_veclist &,
3824 const model_complex_plain_vector &U = md.complex_variable(vl[0]);
3825 return gmm::abs(gmm::vect_hp(matl[0], U, U)) / scalar_type(2);
3829 virtual void asm_complex_tangent_terms(
const model &md,
size_type,
3830 const model::varnamelist &vl,
3831 const model::varnamelist &dl,
3832 const model::mimlist &mims,
3833 model::complex_matlist &matl,
3834 model::complex_veclist &,
3835 model::complex_veclist &,
3837 build_version)
const {
3838 GMM_ASSERT1(matl.size() == 1,
3839 "Generic elliptic brick has one and only one term");
3840 GMM_ASSERT1(mims.size() == 1,
3841 "Generic elliptic brick need one and only one mesh_im");
3842 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
3843 "Wrong number of variables for generic elliptic brick");
3845 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
3846 const mesh &m = mf_u.linked_mesh();
3847 size_type N = m.dim(), Q = mf_u.get_qdim(), s = 1;
3848 const mesh_im &mim = *mims[0];
3849 const model_real_plain_vector *A = 0;
3850 const mesh_fem *mf_a = 0;
3851 mesh_region rg(region);
3852 m.intersect_with_mpi_region(rg);
3855 if (dl.size() > 0) {
3856 A = &(md.real_variable(dl[0]));
3857 mf_a = md.pmesh_fem_of_variable(dl[0]);
3858 s = gmm::vect_size(*A);
3859 if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
3863 GMM_TRACE2(
"Generic elliptic term assembly");
3868 (matl[0], mim, mf_u, *mf_a, *A, rg);
3871 (matl[0], mim, mf_u, *mf_a, *A, rg);
3876 (gmm::real_part(matl[0]), mim, mf_u, rg);
3879 (gmm::real_part(matl[0]), mim, mf_u, rg);
3880 if (A) gmm::scale(matl[0], (*A)[0]);
3882 }
else if (s == N*N) {
3886 (matl[0], mim, mf_u, *mf_a, *A, rg);
3889 (matl[0], mim, mf_u, *mf_a, *A, rg);
3893 (matl[0], mim, mf_u, *A, rg);
3896 (matl[0], mim, mf_u, *A, rg);
3898 }
else if (s == N*N*Q*Q) {
3901 (matl[0], mim, mf_u, *mf_a, *A, rg);
3904 (matl[0], mim, mf_u, *A, rg);
3907 "Bad format generic elliptic brick coefficient");
3910 generic_elliptic_brick() {
3911 set_flags(
"Generic elliptic",
true ,
3919 const std::string &varname,
3922 pbrick pbr = std::make_shared<generic_elliptic_brick>();
3924 tl.push_back(model::term_description(varname, varname,
true));
3925 return md.
add_brick(pbr, model::varnamelist(1, varname),
3926 model::varnamelist(), tl, model::mimlist(1, &mim),
3929 std::string test_varname
3930 =
"Test_" + sup_previous_and_dot_to_varname(varname);
3935 expr =
"Grad_"+varname+
".Grad_"+test_varname;
3937 expr =
"Grad_"+varname+
":Grad_"+test_varname;
3939 "Laplacian",
false);
3944 const std::string &varname,
3945 const std::string &dataname,
3948 pbrick pbr = std::make_shared<generic_elliptic_brick>();
3950 tl.push_back(model::term_description(varname, varname,
true));
3951 return md.
add_brick(pbr, model::varnamelist(1, varname),
3952 model::varnamelist(1, dataname), tl,
3953 model::mimlist(1, &mim), region);
3955 std::string test_varname
3956 =
"Test_" + sup_previous_and_dot_to_varname(varname);
3969 if (qdim_data != 1) {
3970 GMM_ASSERT1(qdim_data == gmm::sqr(dim),
3971 "Wrong data size for generic elliptic brick");
3972 expr =
"((Reshape("+dataname+
",meshdim,meshdim))*Grad_"+varname+
").Grad_"
3975 expr =
"(("+dataname+
")*Grad_"+varname+
").Grad_"+test_varname;
3978 if (qdim_data != 1) {
3979 if (qdim_data == gmm::sqr(dim))
3980 expr =
"((Reshape("+dataname+
",meshdim,meshdim))*Grad_"+varname+
"):Grad_"
3982 else if (qdim_data == gmm::sqr(gmm::sqr(dim)))
3983 expr =
"((Reshape("+dataname+
",meshdim,meshdim,meshdim,meshdim))*Grad_"
3984 +varname+
"):Grad_"+test_varname;
3985 else GMM_ASSERT1(
false,
"Wrong data size for generic elliptic brick");
3987 expr =
"(("+dataname+
")*Grad_"+varname+
"):Grad_"+test_varname;
3991 (md, mim, expr, region,
true,
true,
"Generic elliptic",
true);
3994 "Generic elliptic (nonlinear)");
4006 struct source_term_brick :
public virtual_brick {
4008 void asm_real_tangent_terms(
const model &md,
size_type ,
4009 const model::varnamelist &vl,
4010 const model::varnamelist &dl,
4011 const model::mimlist &mims,
4012 model::real_matlist &,
4013 model::real_veclist &vecl,
4014 model::real_veclist &,
4016 build_version)
const override {
4017 GMM_ASSERT1(vecl.size() == 1,
4018 "Source term brick has one and only one term");
4019 GMM_ASSERT1(mims.size() == 1,
4020 "Source term brick need one and only one mesh_im");
4021 GMM_ASSERT1(vl.size() == 1 && dl.size() > 0 && dl.size() <= 2,
4022 "Wrong number of variables for source term brick");
4024 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4025 const mesh_im &mim = *mims[0];
4026 const model_real_plain_vector &A = md.real_variable(dl[0]);
4027 const mesh_fem *mf_data = md.pmesh_fem_of_variable(dl[0]);
4028 mesh_region rg(region);
4029 mim.linked_mesh().intersect_with_mpi_region(rg);
4032 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4034 GMM_ASSERT1(mf_u.get_qdim() == s,
4035 dl[0] <<
": bad format of source term data. "
4036 "Detected dimension is " << s <<
" should be "
4039 GMM_TRACE2(
"Source term assembly");
4045 if (dl.size() > 1) gmm::add(md.real_variable(dl[1]), vecl[0]);
4048 void real_post_assembly_in_serial(
const model &md,
size_type ib,
4049 const model::varnamelist &,
4050 const model::varnamelist &,
4051 const model::mimlist &,
4052 model::real_matlist &,
4053 model::real_veclist &vecl,
4054 model::real_veclist &,
4056 build_version)
const override
4058 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4062 void asm_complex_tangent_terms(
const model &md,
size_type ,
4063 const model::varnamelist &vl,
4064 const model::varnamelist &dl,
4065 const model::mimlist &mims,
4066 model::complex_matlist &,
4067 model::complex_veclist &vecl,
4068 model::complex_veclist &,
4070 build_version)
const override {
4071 GMM_ASSERT1(vecl.size() == 1,
4072 "Source term brick has one and only one term");
4073 GMM_ASSERT1(mims.size() == 1,
4074 "Source term brick need one and only one mesh_im");
4075 GMM_ASSERT1(vl.size() == 1 && dl.size() > 0 && dl.size() <= 2,
4076 "Wrong number of variables for source term brick");
4078 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4079 const mesh_im &mim = *mims[0];
4080 const model_complex_plain_vector &A = md.complex_variable(dl[0]);
4081 const mesh_fem *mf_data = md.pmesh_fem_of_variable(dl[0]);
4082 mesh_region rg(region);
4083 mim.linked_mesh().intersect_with_mpi_region(rg);
4086 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4088 GMM_ASSERT1(mf_u.get_qdim() == s,
"Bad format of source term data");
4090 GMM_TRACE2(
"Source term assembly");
4096 if (dl.size() > 1)
gmm::add(md.complex_variable(dl[1]), vecl[0]);
4099 void complex_post_assembly_in_serial(
const model &md,
4101 const model::varnamelist &,
4102 const model::varnamelist &,
4103 const model::mimlist &,
4104 model::complex_matlist &,
4105 model::complex_veclist &vecl,
4106 model::complex_veclist &,
4107 size_type, build_version)
const override
4109 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4114 source_term_brick() {
4115 set_flags(
"Source term",
true ,
4125 const std::string &varname,
4126 const std::string &dataexpr,
4128 const std::string &directdataname) {
4130 pbrick pbr = std::make_shared<source_term_brick>();
4132 tl.push_back(model::term_description(varname));
4133 model::varnamelist vdata(1, dataexpr);
4134 if (directdataname.size()) vdata.push_back(directdataname);
4135 return md.
add_brick(pbr, model::varnamelist(1, varname),
4136 vdata, tl, model::mimlist(1, &mim), region);
4138 std::string test_varname
4139 =
"Test_" + sup_previous_and_dot_to_varname(varname);
4144 expr =
"("+dataexpr+
")*"+test_varname;
4146 expr =
"("+dataexpr+
")."+test_varname;
4147 size_type ib = add_source_term_generic_assembly_brick
4148 (md, mim, expr, region,
"Source term", varname, directdataname,
true);
4151 "Source term (nonlinear)");
4152 if (directdataname.size())
4153 add_source_term_generic_assembly_brick
4154 (md, mim,
"", region,
"Source term", varname, directdataname);
4166 struct normal_source_term_brick :
public virtual_brick {
4168 void asm_real_tangent_terms(
const model &md,
size_type ,
4169 const model::varnamelist &vl,
4170 const model::varnamelist &dl,
4171 const model::mimlist &mims,
4172 model::real_matlist &,
4173 model::real_veclist &vecl,
4174 model::real_veclist &,
4176 build_version)
const override {
4177 GMM_ASSERT1(vecl.size() == 1,
4178 "Source term brick has one and only one term");
4179 GMM_ASSERT1(mims.size() == 1,
4180 "Source term brick need one and only one mesh_im");
4181 GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
4182 "Wrong number of variables for source term brick");
4184 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4185 const mesh_im &mim = *mims[0];
4186 const model_real_plain_vector &A = md.real_variable(dl[0]);
4187 const mesh_fem *mf_data = md.pmesh_fem_of_variable(dl[0]);
4188 mesh_region rg(region);
4189 mim.linked_mesh().intersect_with_mpi_region(rg);
4191 size_type s = gmm::vect_size(A), N = mf_u.linked_mesh().dim();
4192 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4194 GMM_ASSERT1(mf_u.get_qdim()*N == s,
4195 dl[0] <<
": bad format of normal source term data. "
4196 "Detected dimension is " << s <<
" should be "
4199 GMM_TRACE2(
"source term assembly");
4206 void real_post_assembly_in_serial(
const model &md,
size_type ib,
4207 const model::varnamelist &,
4208 const model::varnamelist &,
4209 const model::mimlist &,
4210 model::real_matlist &,
4211 model::real_veclist &vecl,
4212 model::real_veclist &,
4214 build_version)
const override {
4215 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4219 virtual void asm_complex_tangent_terms(
const model &md,
size_type ,
4220 const model::varnamelist &vl,
4221 const model::varnamelist &dl,
4222 const model::mimlist &mims,
4223 model::complex_matlist &,
4224 model::complex_veclist &vecl,
4225 model::complex_veclist &,
4227 build_version)
const {
4228 GMM_ASSERT1(vecl.size() == 1,
4229 "Source term brick has one and only one term");
4230 GMM_ASSERT1(mims.size() == 1,
4231 "Source term brick need one and only one mesh_im");
4232 GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
4233 "Wrong number of variables for source term brick");
4235 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4236 const mesh_im &mim = *mims[0];
4237 const model_complex_plain_vector &A = md.complex_variable(dl[0]);
4238 const mesh_fem *mf_data = md.pmesh_fem_of_variable(dl[0]);
4239 mesh_region rg(region);
4240 mim.linked_mesh().intersect_with_mpi_region(rg);
4242 size_type s = gmm::vect_size(A), N = mf_u.linked_mesh().dim();
4243 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4245 GMM_ASSERT1(s == mf_u.get_qdim()*N,
"Bad format of source term data");
4247 GMM_TRACE2(
"source term assembly");
4254 void complex_post_assembly_in_serial(
const model &md,
size_type ib,
4255 const model::varnamelist &,
4256 const model::varnamelist &,
4257 const model::mimlist &,
4258 model::complex_matlist &,
4259 model::complex_veclist &vecl,
4260 model::complex_veclist &,
4262 build_version)
const override {
4263 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4266 normal_source_term_brick() {
4267 set_flags(
"Normal source term",
true ,
4277 const std::string &varname,
4278 const std::string &dataexpr,
4281 pbrick pbr = std::make_shared<normal_source_term_brick>();
4283 tl.push_back(model::term_description(varname));
4284 model::varnamelist vdata(1, dataexpr);
4285 return md.
add_brick(pbr, model::varnamelist(1, varname),
4286 vdata, tl, model::mimlist(1, &mim), region);
4288 std::string test_varname
4289 =
"Test_" + sup_previous_and_dot_to_varname(varname);
4294 expr =
"(("+dataexpr+
").Normal)*"+test_varname;
4296 expr =
"(Reshape("+dataexpr+
",qdim("+varname
4297 +
"),meshdim)*Normal)."+test_varname;
4298 return add_source_term_generic_assembly_brick
4299 (md, mim, expr, region,
"Source term");
4312 struct Dirichlet_condition_brick :
public virtual_brick {
4315 bool normal_component;
4316 const mesh_fem *mf_mult_;
4322 virtual void asm_real_tangent_terms(
const model &md,
size_type ib,
4323 const model::varnamelist &vl,
4324 const model::varnamelist &dl,
4325 const model::mimlist &mims,
4326 model::real_matlist &matl,
4327 model::real_veclist &vecl,
4328 model::real_veclist &,
4330 build_version version)
const {
4331 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
4332 "Dirichlet condition brick has one and only one term");
4333 GMM_ASSERT1(mims.size() == 1,
4334 "Dirichlet condition brick need one and only one mesh_im");
4335 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() <= 3,
4336 "Wrong number of variables for Dirichlet condition brick");
4338 model_real_sparse_matrix& rB = rB_th;
4339 model_real_plain_vector& rV = rV_th;
4341 bool penalized = (vl.size() == 1);
4342 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4343 const mesh_fem &mf_mult = penalized ? (mf_mult_ ? *mf_mult_ : mf_u)
4344 : md.mesh_fem_of_variable(vl[1]);
4345 const mesh_im &mim = *mims[0];
4346 const model_real_plain_vector *A = 0, *COEFF = 0, *H = 0;
4347 const mesh_fem *mf_data = 0, *mf_H = 0;
4348 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
4349 || (penalized && md.is_var_newer_than_brick(dl[0], ib));
4352 COEFF = &(md.real_variable(dl[0]));
4353 GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
4354 "Data for coefficient should be a scalar");
4357 size_type s = 0, ind = (penalized ? 1 : 0);
4358 if (dl.size() > ind) {
4359 A = &(md.real_variable(dl[ind]));
4360 mf_data = md.pmesh_fem_of_variable(dl[ind]);
4361 s = gmm::vect_size(*A);
4362 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4363 size_type ss = ((normal_component) ? 1 : mf_u.get_qdim());
4364 GMM_ASSERT1(s == ss, dl[ind] <<
": bad format of "
4365 "Dirichlet data. Detected dimension is " << s
4366 <<
" should be " << ss);
4369 if (dl.size() > ind + 1) {
4370 GMM_ASSERT1(H_version,
4371 "Wrong number of data for Dirichlet condition brick");
4372 H = &(md.real_variable(dl[ind+1]));
4373 mf_H = md.pmesh_fem_of_variable(dl[ind+1]);
4374 s = gmm::vect_size(*A);
4376 s = s * mf_H->get_qdim() / mf_H->nb_dof();
4377 GMM_ASSERT1(mf_H->get_qdim() == 1,
"Implemented only for mf_H "
4378 "a scalar finite element method");
4380 GMM_ASSERT1(s = gmm::sqr(mf_u.get_qdim()),
4381 dl[ind+1] <<
": bad format of Dirichlet data. "
4382 "Detected dimension is " << s <<
" should be "
4383 <<
size_type(gmm::sqr(mf_u.get_qdim())));
4386 mesh_region rg(region);
4387 mim.linked_mesh().intersect_with_mpi_region(rg);
4389 if (recompute_matrix) {
4390 model_real_sparse_matrix *B = &(matl[0]);
4391 if (penalized && (&mf_mult != &mf_u)) {
4398 GMM_TRACE2(
"Mass term assembly for Dirichlet condition");
4399 if (H_version || normal_component) {
4400 ga_workspace workspace;
4401 gmm::sub_interval Imult(0, mf_mult.nb_dof()), Iu(0, mf_u.nb_dof());
4402 base_vector u(mf_u.nb_dof());
4403 base_vector
mult(mf_mult.nb_dof());
4404 workspace.add_fem_variable(
"u", mf_u, Iu, u);
4405 workspace.add_fem_variable(
"mult", mf_mult, Imult, mult);
4406 auto expression = std::string{};
4408 if (mf_H) workspace.add_fem_constant(
"A", *mf_H, *H);
4409 else workspace.add_fixed_size_constant(
"A", *H);
4410 expression = (mf_u.get_qdim() == 1) ?
4411 "Test_mult . (A . Test2_u)"
4413 "Test_mult. (Reshape(A, qdim(u), qdim(u)) . Test2_u)";
4414 }
else if (normal_component) {
4415 expression =
"Test_mult . (Test2_u . Normal)";
4417 workspace.add_expression(expression, mim, rg);
4418 workspace.set_assembled_matrix(*B);
4419 workspace.assembly(2);
4424 if (penalized && (&mf_mult != &mf_u)) {
4425 GMM_ASSERT1(MPI_IS_MASTER(),
"Sorry, the penalized option "
4426 "filtered by a multiplier space is not parallelized");
4427 gmm::mult(gmm::transposed(rB), rB, matl[0]);
4428 gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
4429 }
else if (penalized) {
4430 gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
4434 if (dl.size() > ind) {
4435 GMM_TRACE2(
"Source term assembly for Dirichlet condition");
4437 if (penalized && (&mf_mult != &mf_u)) {
4451 if (penalized && (&mf_mult != &mf_u)) {
4452 gmm::mult(gmm::transposed(rB), rV, vecl[0]);
4453 gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
4454 rV = model_real_plain_vector();
4455 }
else if (penalized)
4456 gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
4461 void real_post_assembly_in_serial(
const model &md,
size_type ib,
4462 const model::varnamelist &,
4463 const model::varnamelist &,
4464 const model::mimlist &,
4465 model::real_matlist &,
4466 model::real_veclist &vecl,
4467 model::real_veclist &,
4469 build_version)
const override {
4470 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4473 virtual void asm_complex_tangent_terms(
const model &md,
size_type ib,
4474 const model::varnamelist &vl,
4475 const model::varnamelist &dl,
4476 const model::mimlist &mims,
4477 model::complex_matlist &matl,
4478 model::complex_veclist &vecl,
4479 model::complex_veclist &,
4481 build_version version)
const {
4482 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
4483 "Dirichlet condition brick has one and only one term");
4484 GMM_ASSERT1(mims.size() == 1,
4485 "Dirichlet condition brick need one and only one mesh_im");
4486 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() <= 3,
4487 "Wrong number of variables for Dirichlet condition brick");
4489 model_complex_sparse_matrix& cB = cB_th;
4490 model_complex_plain_vector& cV = cV_th;
4492 bool penalized = (vl.size() == 1);
4493 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4494 const mesh_fem &mf_mult = penalized ? (mf_mult_ ? *mf_mult_ : mf_u)
4495 : md.mesh_fem_of_variable(vl[1]);
4496 const mesh_im &mim = *mims[0];
4497 const model_complex_plain_vector *A = 0, *COEFF = 0, *H = 0;
4498 const mesh_fem *mf_data = 0, *mf_H = 0;
4499 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
4500 || (penalized && md.is_var_newer_than_brick(dl[0], ib));
4503 COEFF = &(md.complex_variable(dl[0]));
4504 GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
4505 "Data for coefficient should be a scalar");
4508 size_type s = 0, ind = (penalized ? 1 : 0);
4509 if (dl.size() > ind) {
4510 A = &(md.complex_variable(dl[ind]));
4511 mf_data = md.pmesh_fem_of_variable(dl[ind]);
4512 s = gmm::vect_size(*A);
4513 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4514 size_type ss = s * ((normal_component) ? mf_u.linked_mesh().dim() : 1);
4515 GMM_ASSERT1(mf_u.get_qdim() == ss,
4516 dl[ind] <<
": bad format of Dirichlet data. "
4517 "Detected dimension is " << ss <<
" should be "
4521 if (dl.size() > ind + 1) {
4522 GMM_ASSERT1(H_version,
4523 "Wrong number of data for Dirichlet condition brick");
4524 H = &(md.complex_variable(dl[ind+1]));
4525 mf_H = md.pmesh_fem_of_variable(dl[ind+1]);
4526 s = gmm::vect_size(*A);
4528 s = s * mf_H->get_qdim() / mf_H->nb_dof();
4529 GMM_ASSERT1(mf_H->get_qdim() == 1,
"Implemented only for mf_H "
4530 "a scalar finite element method");
4532 GMM_ASSERT1(s = gmm::sqr(mf_u.get_qdim()),
4533 dl[ind+1] <<
": bad format of Dirichlet data. "
4534 "Detected dimension is " << s <<
" should be "
4535 <<
size_type(gmm::sqr(mf_u.get_qdim())));
4538 mesh_region rg(region);
4539 mim.linked_mesh().intersect_with_mpi_region(rg);
4541 if (recompute_matrix) {
4542 model_complex_sparse_matrix *B = &(matl[0]);
4543 if (penalized && (&mf_mult != &mf_u)) {
4550 GMM_TRACE2(
"Mass term assembly for Dirichlet condition");
4552 if (mf_u.get_qdim() == 1)
4553 asm_real_or_complex_1_param_mat(*B, mim, mf_mult, mf_H, *H, rg,
4554 "(A*Test_u).Test2_u");
4556 asm_real_or_complex_1_param_mat(*B, mim, mf_mult, mf_H, *H, rg,
4557 "(Reshape(A,qdim(u),qdim(u))*Test2_u).Test_u");
4573 else if (normal_component) {
4574 ga_workspace workspace;
4575 gmm::sub_interval Imult(0, mf_mult.nb_dof()), Iu(0, mf_u.nb_dof());
4576 base_vector
mult(mf_mult.nb_dof()), u(mf_u.nb_dof());
4577 workspace.add_fem_variable(
"mult", mf_mult, Imult, mult);
4578 workspace.add_fem_variable(
"u", mf_u, Iu, u);
4579 workspace.add_expression(
"Test_mult.(Test2_u.Normal)", mim, rg);
4580 model_real_sparse_matrix BB(mf_mult.nb_dof(), mf_u.nb_dof());
4581 workspace.set_assembled_matrix(BB);
4582 workspace.assembly(2);
4598 if (penalized && (&mf_mult != &mf_u)) {
4599 gmm::mult(gmm::transposed(cB), cB, matl[0]);
4600 gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
4601 }
else if (penalized) {
4602 gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
4606 if (dl.size() > ind) {
4607 GMM_TRACE2(
"Source term assembly for Dirichlet condition");
4609 if (penalized && (&mf_mult != &mf_u)) {
4623 if (penalized && (&mf_mult != &mf_u)) {
4624 gmm::mult(gmm::transposed(cB), cV, vecl[0]);
4625 gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
4626 cV = model_complex_plain_vector();
4627 }
else if (penalized)
4628 gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
4632 void complex_post_assembly_in_serial(
const model &md,
size_type ib,
4633 const model::varnamelist &,
4634 const model::varnamelist &,
4635 const model::mimlist &,
4636 model::complex_matlist &,
4637 model::complex_veclist &vecl,
4638 model::complex_veclist &,
4640 build_version)
const override {
4641 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4645 virtual std::string declare_volume_assembly_string
4646 (
const model &,
size_type,
const model::varnamelist &,
4647 const model::varnamelist &)
const {
4648 return std::string();
4651 Dirichlet_condition_brick(
bool penalized,
bool H_version_,
4652 bool normal_component_,
4653 const mesh_fem *mf_mult__ = 0) {
4654 mf_mult_ = mf_mult__;
4655 H_version = H_version_;
4656 normal_component = normal_component_;
4657 GMM_ASSERT1(!(H_version && normal_component),
"Bad Dirichlet version");
4658 set_flags(penalized ?
"Dirichlet with penalization brick"
4659 :
"Dirichlet with multipliers brick",
4668 (
model &md,
const mesh_im &mim,
const std::string &varname,
4669 const std::string &multname,
size_type region,
4670 const std::string &dataname) {
4671 pbrick pbr = std::make_shared<Dirichlet_condition_brick>(
false,
false,
false);
4673 tl.push_back(model::term_description(multname, varname,
true));
4674 model::varnamelist vl(1, varname);
4675 vl.push_back(multname);
4676 model::varnamelist dl;
4677 if (dataname.size()) dl.push_back(dataname);
4678 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4682 (
model &md,
const mesh_im &mim,
const std::string &varname,
4684 const std::string &dataname) {
4685 std::string multname = md.
new_name(
"mult_on_" + varname);
4688 (md, mim, varname, multname, region, dataname);
4692 (
model &md,
const mesh_im &mim,
const std::string &varname,
4694 const std::string &dataname) {
4699 (md, mim, varname, mf_mult, region, dataname);
4707 (
model &md,
const mesh_im &mim,
const std::string &varname,
4708 scalar_type penalisation_coeff,
size_type region,
4709 const std::string &dataname,
const mesh_fem *mf_mult) {
4710 std::string coeffname = md.
new_name(
"penalization_on_" + varname);
4716 pbrick pbr = std::make_shared<Dirichlet_condition_brick>
4717 (
true,
false,
false, mf_mult);
4719 tl.push_back(model::term_description(varname, varname,
true));
4720 model::varnamelist vl(1, varname);
4721 model::varnamelist dl(1, coeffname);
4722 if (dataname.size()) dl.push_back(dataname);
4723 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4727 (
model &md,
const mesh_im &mim,
const std::string &varname,
4728 const std::string &multname,
size_type region,
4729 const std::string &dataname) {
4730 pbrick pbr = std::make_shared<Dirichlet_condition_brick>(
false,
false,
true);
4732 tl.push_back(model::term_description(multname, varname,
true));
4733 model::varnamelist vl(1, varname);
4734 vl.push_back(multname);
4735 model::varnamelist dl;
4736 if (dataname.size()) dl.push_back(dataname);
4737 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4741 (
model &md,
const mesh_im &mim,
const std::string &varname,
4743 const std::string &dataname) {
4744 std::string multname = md.
new_name(
"mult_on_" + varname);
4747 (md, mim, varname, multname, region, dataname);
4751 (
model &md,
const mesh_im &mim,
const std::string &varname,
4753 const std::string &dataname) {
4757 (md, mim, varname, mf_mult, region, dataname);
4761 (
model &md,
const mesh_im &mim,
const std::string &varname,
4762 scalar_type penalisation_coeff,
size_type region,
4763 const std::string &dataname,
const mesh_fem *mf_mult) {
4764 std::string coeffname = md.
new_name(
"penalization_on_" + varname);
4770 pbrick pbr = std::make_shared<Dirichlet_condition_brick>
4771 (
true,
false,
true, mf_mult);
4773 tl.push_back(model::term_description(varname, varname,
true));
4774 model::varnamelist vl(1, varname);
4775 model::varnamelist dl(1, coeffname);
4776 if (dataname.size()) dl.push_back(dataname);
4777 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4782 (
model &md,
const mesh_im &mim,
const std::string &varname,
4783 const std::string &multname,
size_type region,
4784 const std::string &dataname,
const std::string &Hname) {
4785 pbrick pbr = std::make_shared<Dirichlet_condition_brick>(
false,
true,
false);
4787 tl.push_back(model::term_description(multname, varname,
true));
4788 model::varnamelist vl(1, varname);
4789 vl.push_back(multname);
4790 model::varnamelist dl;
4791 dl.push_back(dataname);
4792 dl.push_back(Hname);
4793 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4797 (
model &md,
const mesh_im &mim,
const std::string &varname,
4799 const std::string &dataname,
const std::string &Hname) {
4800 std::string multname = md.
new_name(
"mult_on_" + varname);
4803 (md, mim, varname, multname, region, dataname, Hname);
4807 (
model &md,
const mesh_im &mim,
const std::string &varname,
4809 const std::string &dataname,
const std::string &Hname) {
4814 (md, mim, varname, mf_mult, region, dataname, Hname);
4818 (
model &md,
const mesh_im &mim,
const std::string &varname,
4819 scalar_type penalisation_coeff,
size_type region,
4820 const std::string &dataname,
const std::string &Hname,
4822 std::string coeffname = md.
new_name(
"penalization_on_" + varname);
4828 pbrick pbr = std::make_shared<Dirichlet_condition_brick>
4829 (
true,
true,
false, mf_mult);
4831 tl.push_back(model::term_description(varname, varname,
true));
4832 model::varnamelist vl(1, varname);
4833 model::varnamelist dl(1, coeffname);
4834 dl.push_back(dataname);
4835 dl.push_back(Hname);
4836 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4840 scalar_type penalisation_coeff) {
4844 GMM_ASSERT1(gmm::vect_size(d)==1,
4845 "Wrong coefficient size, may be not a Dirichlet brick "
4846 "with penalization");
4847 d[0] = penalisation_coeff;
4851 GMM_ASSERT1(gmm::vect_size(d)==1,
4852 "Wrong coefficient size, may be not a Dirichlet brick "
4853 "with penalization");
4854 d[0] = penalisation_coeff;
4864 struct simplification_Dirichlet_condition_brick :
public virtual_brick {
4866 virtual void real_pre_assembly_in_serial(
const model &md,
size_type ,
4867 const model::varnamelist &vl,
4868 const model::varnamelist &dl,
4869 const model::mimlist &mims,
4870 model::real_matlist &matl,
4871 model::real_veclist &vecl,
4872 model::real_veclist &,
4874 build_version )
const {
4875 if (MPI_IS_MASTER()) {
4877 GMM_ASSERT1(vecl.size() == 0 && matl.size() == 0,
4878 "Dirichlet condition brick by simplification has no term");
4879 GMM_ASSERT1(mims.size() == 0,
4880 "Dirichlet condition brick by simplification need no "
4882 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
4883 "Wrong number of variables for Dirichlet condition brick "
4884 "by simplification");
4886 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4887 const model_real_plain_vector *A = 0;
4888 const mesh_fem *mf_data = 0;
4891 if (dl.size() == 1) {
4892 A = &(md.real_variable(dl[0]));
4893 mf_data = md.pmesh_fem_of_variable(dl[0]);
4896 GMM_ASSERT1(mf_data == &mf_u,
"Sorry, for this brick, the data has"
4897 " to be defined on the same f.e.m. as the unknown");
4899 s = gmm::vect_size(*A);
4900 GMM_ASSERT1(mf_u.get_qdim() == s,
": bad format of "
4901 "Dirichlet data. Detected dimension is " << s
4902 <<
" should be " <<
size_type(mf_u.get_qdim()));
4906 mesh_region rg(region);
4909 if (mf_u.get_qdim() > 1 || (!mf_data && A)) {
4910 for (mr_visitor i(rg, mf_u.linked_mesh()); !i.finished(); ++i) {
4911 pfem pf = mf_u.fem_of_element(i.cv());
4913 GMM_ASSERT1(pf->target_dim() == 1,
4914 "Intrinsically vectorial fems are not allowed");
4915 GMM_ASSERT1(mf_data || pf->is_lagrange(),
"Constant Dirichlet "
4916 "data allowed for lagrange fems only");
4921 dal::bit_vector dofs = mf_u.dof_on_region(rg);
4923 if (A && !mf_data) {
4924 GMM_ASSERT1(dofs.card() % s == 0,
"Problem with dof vectorization");
4927 for (dal::bv_visitor i(dofs); !i.finished(); ++i) {
4929 if (A) val = (mf_data ? (*A)[i] : (*A)[i%s]);
4930 md.add_real_dof_constraint(vl[0], i, val);
4935 virtual void complex_pre_assembly_in_serial(
const model &md,
size_type ,
4936 const model::varnamelist &vl,
4937 const model::varnamelist &dl,
4938 const model::mimlist &mims,
4939 model::complex_matlist &matl,
4940 model::complex_veclist &vecl,
4941 model::complex_veclist &,
4943 build_version )
const {
4944 if (MPI_IS_MASTER()) {
4945 GMM_ASSERT1(vecl.size() == 0 && matl.size() == 0,
4946 "Dirichlet condition brick by simplification has no term");
4947 GMM_ASSERT1(mims.size() == 0,
4948 "Dirichlet condition brick by simplification need no "
4950 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
4951 "Wrong number of variables for Dirichlet condition brick "
4952 "by simplification");
4954 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4955 const model_complex_plain_vector *A = 0;
4956 const mesh_fem *mf_data = 0;
4959 if (dl.size() == 1) {
4960 A = &(md.complex_variable(dl[0]));
4961 mf_data = md.pmesh_fem_of_variable(dl[0]);
4964 GMM_ASSERT1(mf_data == &mf_u,
"Sorry, for this brick, the data has"
4965 " to be define on the same f.e.m. than the unknown");
4967 s = gmm::vect_size(*A);
4968 GMM_ASSERT1(mf_u.get_qdim() == s,
": bad format of "
4969 "Dirichlet data. Detected dimension is " << s
4970 <<
" should be " <<
size_type(mf_u.get_qdim()));
4974 mesh_region rg(region);
4977 if (mf_u.get_qdim() > 1 || (!mf_data && A)) {
4978 for (mr_visitor i(rg, mf_u.linked_mesh()); !i.finished(); ++i) {
4979 pfem pf = mf_u.fem_of_element(i.cv());
4981 GMM_ASSERT1(pf->target_dim() == 1,
4982 "Intrinsically vectorial fems are not allowed");
4983 GMM_ASSERT1(mf_data || pf->is_lagrange(),
"Constant Dirichlet "
4984 "data allowed for lagrange fems only");
4989 dal::bit_vector dofs = mf_u.dof_on_region(rg);
4991 if (A && !mf_data) {
4992 GMM_ASSERT1(dofs.card() % s == 0,
"Problem with dof vectorization");
4995 for (dal::bv_visitor i(dofs); !i.finished(); ++i) {
4996 complex_type val(0);
4997 if (A) val = (mf_data ? (*A)[i] : (*A)[i%s]);
4998 md.add_complex_dof_constraint(vl[0], i, val);
5004 virtual std::string declare_volume_assembly_string
5005 (
const model &,
size_type,
const model::varnamelist &,
5006 const model::varnamelist &)
const {
5007 return std::string();
5010 simplification_Dirichlet_condition_brick() {
5011 set_flags(
"Dirichlet with simplification brick",
5020 (
model &md,
const std::string &varname,
5021 size_type region,
const std::string &dataname) {
5022 pbrick pbr = std::make_shared<simplification_Dirichlet_condition_brick>();
5024 model::varnamelist vl(1, varname);
5025 model::varnamelist dl;
5026 if (dataname.size()) dl.push_back(dataname);
5027 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(), region);
5037 (
model &md,
const mesh_im &mim,
const std::string &varname,
5038 const std::string &Neumannterm,
5039 const std::string &datagamma0,
size_type region, scalar_type theta_,
5040 const std::string &datag) {
5041 std::string theta = std::to_string(theta_);
5043 ga_workspace workspace(md, ga_workspace::inherit::ALL);
5044 size_type order = workspace.add_expression(Neumannterm, mim, region, 1);
5045 GMM_ASSERT1(order == 0,
"Wrong expression of the Neumann term");
5046 bool is_lin = workspace.is_linear(1);
5048 std::string condition =
"("+varname + (datag.size() ?
"-("+datag+
"))":
")");
5049 std::string gamma =
"(("+datagamma0+
")*element_size)";
5050 std::string r =
"(1/"+gamma+
")";
5051 std::string expr =
"("+r+
"*"+condition+
"-("+Neumannterm+
")).Test_"+varname;
5052 if (theta_ != scalar_type(0)) {
5053 std::string derivative_Neumann = workspace.extract_order1_term(varname);
5054 if (derivative_Neumann.size())
5055 expr+=
"-"+theta+
"*"+condition+
".("+derivative_Neumann+
")";
5063 "Dirichlet condition with Nitsche's method");
5066 "Dirichlet condition with Nitsche's method");
5071 (
model &md,
const mesh_im &mim,
const std::string &varname,
5072 const std::string &Neumannterm,
5073 const std::string &datagamma0,
size_type region, scalar_type theta_,
5074 const std::string &datag) {
5075 std::string theta = std::to_string(theta_);
5077 ga_workspace workspace(md, ga_workspace::inherit::ALL);
5078 size_type order = workspace.add_expression(Neumannterm, mim, region, 1);
5079 GMM_ASSERT1(order == 0,
"Wrong expression of the Neumann term");
5080 bool is_lin = workspace.is_linear(1);
5082 std::string condition =
"("+varname+
".Normal"
5083 + (datag.size() ?
"-("+datag+
"))":
")");
5084 std::string gamma =
"(("+datagamma0+
")*element_size)";
5085 std::string r =
"(1/"+gamma+
")";
5086 std::string expr =
"("+r+
"*"+condition+
"-Normal.("+Neumannterm
5087 +
"))*(Normal.Test_"+varname+
")";
5088 if (theta_ != scalar_type(0)) {
5089 std::string derivative_Neumann = workspace.extract_order1_term(varname);
5090 if (derivative_Neumann.size())
5091 expr+=
"-"+theta+
"*"+condition+
"*Normal.("
5092 +derivative_Neumann+
")";
5096 "Dirichlet condition with Nitsche's method");
5099 "Dirichlet condition with Nitsche's method");
5104 (
model &md,
const mesh_im &mim,
const std::string &varname,
5105 const std::string &Neumannterm,
5106 const std::string &datagamma0,
size_type region, scalar_type theta_,
5107 const std::string &datag,
const std::string &dataH) {
5108 std::string theta = std::to_string(theta_);
5110 ga_workspace workspace(md, ga_workspace::inherit::ALL);
5111 size_type order = workspace.add_expression(Neumannterm, mim, region, 1);
5112 GMM_ASSERT1(order == 0,
"Wrong expression of the Neumann term");
5113 bool is_lin = workspace.is_linear(1);
5115 std::string condition =
"(("+dataH+
")*"+varname
5116 + (datag.size() ?
"-("+datag+
"))":
")");
5117 std::string gamma =
"(("+datagamma0+
")*element_size)";
5118 std::string r =
"(1/"+gamma+
")";
5119 std::string expr =
"("+r+
"*"+condition+
"-("+dataH+
")*("+Neumannterm
5120 +
"))*(("+dataH+
")*Test_"+varname+
")";
5121 if (theta_ != scalar_type(0)) {
5122 std::string derivative_Neumann = workspace.extract_order1_term(varname);
5123 if (derivative_Neumann.size())
5124 expr+=
"-"+theta+
"*"+condition+
"*(("+dataH+
")*("
5125 +derivative_Neumann+
"))";
5129 "Dirichlet condition with Nitsche's method");
5132 "Dirichlet condition with Nitsche's method");
5144 struct pointwise_constraints_brick :
public virtual_brick {
5146 mutable gmm::row_matrix<model_real_sparse_vector> rB;
5147 mutable gmm::row_matrix<model_complex_sparse_vector> cB;
5149 virtual void real_pre_assembly_in_serial(
const model &md,
size_type ib,
5150 const model::varnamelist &vl,
5151 const model::varnamelist &dl,
5152 const model::mimlist &mims,
5153 model::real_matlist &matl,
5154 model::real_veclist &vecl,
5155 model::real_veclist &,
5157 build_version version)
const {
5158 if (MPI_IS_MASTER()) {
5160 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5161 "Pointwize constraints brick has only one term");
5162 GMM_ASSERT1(mims.size() == 0,
5163 "Pointwize constraints brick does not need a mesh_im");
5164 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2,
5165 "Wrong number of variables for pointwize constraints brick");
5166 bool penalized = (vl.size() == 1);
5167 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5168 dim_type N = mf_u.linked_mesh().dim(), Q = mf_u.get_qdim(), ind_pt = 0;
5170 GMM_ASSERT1(dl.size() == dlsize || dl.size() == dlsize+1,
5171 "Wrong number of data for pointwize constraints brick");
5174 const model_real_plain_vector *COEFF = 0;
5176 COEFF = &(md.real_variable(dl[0]));
5178 GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
5179 "Data for coefficient should be a scalar");
5182 const model_real_plain_vector &PT = md.real_variable(dl[ind_pt]);
5183 size_type nb_co = gmm::vect_size(PT) / N;
5185 dim_type ind_unitv = dim_type((Q > 1) ? ind_pt+1 : 0);
5186 const model_real_plain_vector &unitv =md.real_variable(dl[ind_unitv]);
5187 GMM_ASSERT1((!ind_unitv || gmm::vect_size(unitv) == nb_co * Q),
5188 "Wrong size for vector of unit vectors");
5190 dim_type ind_rhs = dim_type((Q > 1) ? ind_pt+2 : ind_pt+1);
5191 if (dl.size() <
size_type(ind_rhs + 1)) ind_rhs = 0;
5192 const model_real_plain_vector &rhs = md.real_variable(dl[ind_rhs]);
5193 GMM_ASSERT1((!ind_rhs || gmm::vect_size(rhs) == nb_co),
5194 "Wrong size for vector of rhs");
5196 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
5197 || (penalized && (md.is_var_newer_than_brick(dl[ind_pt], ib)
5198 || md.is_var_newer_than_brick(dl[ind_unitv], ib)
5199 || md.is_var_newer_than_brick(dl[ind_rhs], ib)));
5201 if (recompute_matrix) {
5202 gmm::row_matrix<model_real_sparse_vector> BB(nb_co*Q, mf_u.nb_dof());
5205 dal::bit_vector dof_untouched;
5206 getfem::mesh_trans_inv mti(mf_u.linked_mesh());
5209 gmm::copy(gmm::sub_vector(PT, gmm::sub_interval(i*N, N)), pt);
5212 gmm::row_matrix<model_real_sparse_vector> &BBB = ((Q > 1) ? BB : rB);
5213 model_real_plain_vector vv;
5214 interpolation(mf_u, mti, vv, vv, BBB, 1, 1, &dof_untouched);
5215 GMM_ASSERT1(dof_untouched.card() == 0,
5216 "Pointwize constraints : some of the points are outside "
5217 "the mesh: " << dof_untouched);
5222 gmm::add(gmm::scaled(gmm::mat_row(BB, i*Q+q), unitv[i*Q+q]),
5223 gmm::mat_row(rB, i));
5226 gmm::mult(gmm::transposed(rB), rB, matl[0]);
5227 gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
5234 gmm::mult(gmm::transposed(rB), rhs, vecl[0]);
5235 gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
5243 virtual void complex_pre_assembly_in_serial(
const model &md,
size_type ib,
5244 const model::varnamelist &vl,
5245 const model::varnamelist &dl,
5246 const model::mimlist &mims,
5247 model::complex_matlist &matl,
5248 model::complex_veclist &vecl,
5249 model::complex_veclist &,
5251 build_version version)
const {
5252 if (MPI_IS_MASTER()) {
5253 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5254 "Pointwize constraints brick only one term");
5255 GMM_ASSERT1(mims.size() == 0,
5256 "Pointwize constraints brick does not need a mesh_im");
5257 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2,
5258 "Wrong number of variables for pointwize constraints brick");
5259 bool penalized = (vl.size() == 1);
5260 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5261 dim_type N = mf_u.linked_mesh().dim(), Q = mf_u.get_qdim(), ind_pt = 0;
5263 GMM_ASSERT1(dl.size() == dlsize || dl.size() == dlsize+1,
5264 "Wrong number of data for pointwize constraints brick");
5267 const model_complex_plain_vector *COEFF = 0;
5269 COEFF = &(md.complex_variable(dl[0]));
5271 GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
5272 "Data for coefficient should be a scalar");
5275 const model_complex_plain_vector &PT = md.complex_variable(dl[ind_pt]);
5276 size_type nb_co = gmm::vect_size(PT) / N;
5278 dim_type ind_unitv = dim_type((Q > 1) ? ind_pt+1 : 0);
5279 const model_complex_plain_vector &unitv
5280 = md.complex_variable(dl[ind_unitv]);
5281 GMM_ASSERT1((!ind_unitv || gmm::vect_size(unitv) == nb_co * Q),
5282 "Wrong size for vector of unit vectors");
5284 dim_type ind_rhs = dim_type((Q > 1) ? ind_pt+2 : ind_pt+1);
5285 if (dl.size() <
size_type(ind_rhs + 1)) ind_rhs = 0;
5286 const model_complex_plain_vector &rhs
5287 = md.complex_variable(dl[ind_rhs]);
5288 GMM_ASSERT1((!ind_rhs || gmm::vect_size(rhs) == nb_co),
5289 "Wrong size for vector of rhs");
5291 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
5292 || (penalized && (md.is_var_newer_than_brick(dl[ind_pt], ib)
5293 || md.is_var_newer_than_brick(dl[ind_unitv], ib)
5294 || md.is_var_newer_than_brick(dl[ind_rhs], ib)));
5296 if (recompute_matrix) {
5297 gmm::row_matrix<model_complex_sparse_vector>
5298 BB(nb_co*Q,mf_u.nb_dof());
5300 dal::bit_vector dof_untouched;
5301 getfem::mesh_trans_inv mti(mf_u.linked_mesh());
5305 (gmm::sub_vector(PT, gmm::sub_interval(i*N, N))), pt);
5308 gmm::row_matrix<model_complex_sparse_vector> &BBB = ((Q > 1) ? BB :cB);
5309 model_complex_plain_vector vv;
5310 interpolation(mf_u, mti, vv, vv, BBB, 1, 1, &dof_untouched);
5311 GMM_ASSERT1(dof_untouched.card() == 0,
5312 "Pointwize constraints : some of the points are outside "
5313 "the mesh: " << dof_untouched);
5318 gmm::add(gmm::scaled(gmm::mat_row(BB, i*Q+q), unitv[i*Q+q]),
5319 gmm::mat_row(cB, i));
5323 gmm::mult(gmm::transposed(cB), cB, matl[0]);
5324 gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
5332 gmm::mult(gmm::transposed(cB), rhs, vecl[0]);
5333 gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
5341 virtual std::string declare_volume_assembly_string
5342 (
const model &,
size_type,
const model::varnamelist &,
5343 const model::varnamelist &)
const {
5344 return std::string();
5347 pointwise_constraints_brick(
bool penalized) {
5348 set_flags(penalized ?
"Pointwise cosntraints with penalization brick"
5349 :
"Pointwise cosntraints with multipliers brick",
5359 (
model &md,
const std::string &varname,
5360 scalar_type penalisation_coeff,
const std::string &dataname_pt,
5361 const std::string &dataname_unitv,
const std::string &dataname_val) {
5362 std::string coeffname = md.
new_name(
"penalization_on_" + varname);
5368 pbrick pbr = std::make_shared<pointwise_constraints_brick>(
true);
5370 tl.push_back(model::term_description(varname, varname,
true));
5371 model::varnamelist vl(1, varname);
5372 model::varnamelist dl(1, coeffname);
5373 dl.push_back(dataname_pt);
5375 if (mf_u.
get_qdim() > 1) dl.push_back(dataname_unitv);
5376 if (dataname_val.size() > 0) dl.push_back(dataname_val);
5381 (
model &md,
const std::string &varname,
5382 const std::string &multname,
const std::string &dataname_pt,
5383 const std::string &dataname_unitv,
const std::string &dataname_val) {
5384 pbrick pbr = std::make_shared<pointwise_constraints_brick>(
false);
5386 tl.push_back(model::term_description(multname, varname,
true));
5387 model::varnamelist vl(1, varname);
5388 vl.push_back(multname);
5389 model::varnamelist dl(1, dataname_pt);
5391 if (mf_u.
get_qdim() > 1) dl.push_back(dataname_unitv);
5392 if (dataname_val.size() > 0) dl.push_back(dataname_val);
5397 (
model &md,
const std::string &varname,
const std::string &dataname_pt,
5398 const std::string &dataname_unitv,
const std::string &dataname_val) {
5399 std::string multname = md.
new_name(
"mult_on_" + varname);
5407 (md, varname, multname, dataname_pt, dataname_unitv, dataname_val);
5417 struct Helmholtz_brick :
public virtual_brick {
5419 virtual void asm_real_tangent_terms(
const model &md,
size_type,
5420 const model::varnamelist &vl,
5421 const model::varnamelist &dl,
5422 const model::mimlist &mims,
5423 model::real_matlist &matl,
5424 model::real_veclist &,
5425 model::real_veclist &,
5427 build_version)
const {
5428 GMM_ASSERT1(matl.size() == 1,
5429 "Helmholtz brick has one and only one term");
5430 GMM_ASSERT1(mims.size() == 1,
5431 "Helmholtz brick need one and only one mesh_im");
5432 GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
5433 "Wrong number of variables for Helmholtz brick");
5435 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5436 const mesh &m = mf_u.linked_mesh();
5438 GMM_ASSERT1(Q == 1,
"Helmholtz brick is only for scalar field, sorry.");
5439 const mesh_im &mim = *mims[0];
5440 const mesh_fem *mf_a = 0;
5441 mesh_region rg(region);
5442 m.intersect_with_mpi_region(rg);
5443 const model_real_plain_vector *A = &(md.real_variable(dl[0]));
5444 mf_a = md.pmesh_fem_of_variable(dl[0]);
5445 s = gmm::vect_size(*A);
5446 if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
5449 GMM_TRACE2(
"Stiffness matrix assembly for Helmholtz problem");
5450 gmm::clear(matl[0]);
5451 model_real_plain_vector A2(gmm::vect_size(*A));
5452 for (
size_type i=0; i < gmm::vect_size(*A); ++i)
5453 A2[i] = gmm::sqr((*A)[i]);
5459 GMM_ASSERT1(
false,
"Bad format Helmholtz brick coefficient");
5462 virtual void asm_complex_tangent_terms(
const model &md,
size_type,
5463 const model::varnamelist &vl,
5464 const model::varnamelist &dl,
5465 const model::mimlist &mims,
5466 model::complex_matlist &matl,
5467 model::complex_veclist &,
5468 model::complex_veclist &,
5470 build_version)
const {
5471 GMM_ASSERT1(matl.size() == 1,
5472 "Helmholtz brick has one and only one term");
5473 GMM_ASSERT1(mims.size() == 1,
5474 "Helmholtz brick need one and only one mesh_im");
5475 GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
5476 "Wrong number of variables for Helmholtz brick");
5478 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5479 const mesh &m = mf_u.linked_mesh();
5481 GMM_ASSERT1(Q == 1,
"Helmholtz brick is only for scalar field, sorry.");
5482 const mesh_im &mim = *mims[0];
5483 const mesh_fem *mf_a = 0;
5484 mesh_region rg(region);
5485 m.intersect_with_mpi_region(rg);
5486 const model_complex_plain_vector *A = &(md.complex_variable(dl[0]));
5487 mf_a = md.pmesh_fem_of_variable(dl[0]);
5488 s = gmm::vect_size(*A);
5489 if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
5492 GMM_TRACE2(
"Stiffness matrix assembly for Helmholtz problem");
5494 model_complex_plain_vector A2(gmm::vect_size(*A));
5495 for (
size_type i=0; i < gmm::vect_size(*A); ++i)
5496 A2[i] = gmm::sqr((*A)[i]);
5502 GMM_ASSERT1(
false,
"Bad format Helmholtz brick coefficient");
5506 set_flags(
"Helmholtz",
true ,
5514 const std::string &varname,
5515 const std::string &dataexpr,
5518 pbrick pbr = std::make_shared<Helmholtz_brick>();
5520 tl.push_back(model::term_description(varname, varname,
true));
5521 return md.
add_brick(pbr, model::varnamelist(1, varname),
5522 model::varnamelist(1, dataexpr), tl,
5523 model::mimlist(1, &mim), region);
5525 std::string test_varname
5526 =
"Test_" + sup_previous_and_dot_to_varname(varname);
5527 std::string expr =
"Grad_"+varname+
".Grad_"+test_varname
5528 +
" + sqr("+dataexpr+
")*"+varname+
"*"+test_varname;
5534 "Helmholtz (nonlinear)");
5547 struct Fourier_Robin_brick :
public virtual_brick {
5549 virtual void asm_real_tangent_terms(
const model &md,
size_type,
5550 const model::varnamelist &vl,
5551 const model::varnamelist &dl,
5552 const model::mimlist &mims,
5553 model::real_matlist &matl,
5554 model::real_veclist &,
5555 model::real_veclist &,
5557 build_version)
const {
5558 GMM_ASSERT1(matl.size() == 1,
5559 "Fourier-Robin brick has one and only one term");
5560 GMM_ASSERT1(mims.size() == 1,
5561 "Fourier-Robin brick need one and only one mesh_im");
5562 GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
5563 "Wrong number of variables for Fourier-Robin brick");
5565 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5566 const mesh &m = mf_u.linked_mesh();
5568 const mesh_im &mim = *mims[0];
5569 const mesh_fem *mf_a = 0;
5570 mesh_region rg(region);
5571 m.intersect_with_mpi_region(rg);
5572 const model_real_plain_vector *A = &(md.real_variable(dl[0]));
5573 mf_a = md.pmesh_fem_of_variable(dl[0]);
5574 s = gmm::vect_size(*A);
5575 if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
5576 GMM_ASSERT1(s == Q*Q,
5577 "Bad format Fourier-Robin brick coefficient");
5579 GMM_TRACE2(
"Fourier-Robin term assembly");
5580 gmm::clear(matl[0]);
5584 asm_homogeneous_qu_term(matl[0], mim, mf_u, *A, rg);
5587 virtual void asm_complex_tangent_terms(
const model &md,
size_type,
5588 const model::varnamelist &vl,
5589 const model::varnamelist &dl,
5590 const model::mimlist &mims,
5591 model::complex_matlist &matl,
5592 model::complex_veclist &,
5593 model::complex_veclist &,
5595 build_version)
const {
5596 GMM_ASSERT1(matl.size() == 1,
5597 "Fourier-Robin brick has one and only one term");
5598 GMM_ASSERT1(mims.size() == 1,
5599 "Fourier-Robin brick need one and only one mesh_im");
5600 GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
5601 "Wrong number of variables for Fourier-Robin brick");
5603 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5604 const mesh &m = mf_u.linked_mesh();
5606 const mesh_im &mim = *mims[0];
5607 const mesh_fem *mf_a = 0;
5608 mesh_region rg(region);
5609 m.intersect_with_mpi_region(rg);
5610 const model_complex_plain_vector *A = &(md.complex_variable(dl[0]));
5611 mf_a = md.pmesh_fem_of_variable(dl[0]);
5612 s = gmm::vect_size(*A);
5613 if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
5614 GMM_ASSERT1(s == Q*Q,
5615 "Bad format Fourier-Robin brick coefficient");
5617 GMM_TRACE2(
"Fourier-Robin term assembly");
5622 asm_homogeneous_qu_term(matl[0], mim, mf_u, *A, rg);
5625 Fourier_Robin_brick() {
5626 set_flags(
"Fourier Robin condition",
true ,
5635 const std::string &varname,
5636 const std::string &dataexpr,
5639 pbrick pbr = std::make_shared<Fourier_Robin_brick>();
5641 tl.push_back(model::term_description(varname, varname,
true));
5642 return md.
add_brick(pbr, model::varnamelist(1, varname),
5643 model::varnamelist(1, dataexpr), tl,
5644 model::mimlist(1, &mim), region);
5646 std::string test_varname
5647 =
"Test_" + sup_previous_and_dot_to_varname(varname);
5648 std::string expr =
"(("+dataexpr+
")*"+varname+
")."+test_varname;
5650 "Fourier-Robin",
true);
5653 "Fourier-Robin (nonlinear)");
5664 struct have_private_data_brick :
public virtual_brick {
5666 model_real_sparse_matrix rB;
5667 model_complex_sparse_matrix cB;
5668 model_real_plain_vector rL;
5669 model_complex_plain_vector cL;
5673 struct constraint_brick :
public have_private_data_brick {
5675 virtual void real_pre_assembly_in_serial(
const model &md,
size_type,
5676 const model::varnamelist &vl,
5677 const model::varnamelist &dl,
5678 const model::mimlist &mims,
5679 model::real_matlist &matl,
5680 model::real_veclist &vecl,
5681 model::real_veclist &,
5683 if (MPI_IS_MASTER()) {
5685 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5686 "Constraint brick has one and only one term");
5687 GMM_ASSERT1(mims.size() == 0,
5688 "Constraint brick need no mesh_im");
5689 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() <= 1,
5690 "Wrong number of variables for constraint brick");
5692 bool penalized = (vl.size() == 1);
5693 const model_real_plain_vector *COEFF = 0;
5695 bool has_data = (nameL.compare(
"") != 0);
5697 GMM_ASSERT1(nameL.compare(dl.back()) == 0 &&
5698 md.variable_exists(nameL) && md.is_data(nameL),
5700 const model_real_plain_vector &
5701 rrL = has_data ? md.real_variable(nameL) : rL;
5704 COEFF = &(md.real_variable(dl[0]));
5705 GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
5706 "Data for coefficient should be a scalar");
5709 gmm::scaled(rrL, gmm::abs((*COEFF)[0])), vecl[0]);
5711 gmm::scaled(rB, gmm::abs((*COEFF)[0])), matl[0]);
5719 virtual void complex_pre_assembly_in_serial(
const model &md,
size_type,
5720 const model::varnamelist &vl,
5721 const model::varnamelist &dl,
5722 const model::mimlist &mims,
5723 model::complex_matlist &matl,
5724 model::complex_veclist &vecl,
5725 model::complex_veclist &,
5727 build_version)
const {
5728 if (MPI_IS_MASTER()) {
5730 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5731 "Constraint brick has one and only one term");
5732 GMM_ASSERT1(mims.size() == 0,
5733 "Constraint brick need no mesh_im");
5734 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() <= 1,
5735 "Wrong number of variables for constraint brick");
5737 bool penalized = (vl.size() == 1);
5738 const model_complex_plain_vector *COEFF = 0;
5740 bool has_data = (nameL.compare(
"") != 0);
5742 GMM_ASSERT1(nameL.compare(dl.back()) == 0 &&
5743 md.variable_exists(nameL) && md.is_data(nameL),
5745 const model_complex_plain_vector &
5746 ccL = has_data ? md.complex_variable(nameL) : cL;
5749 COEFF = &(md.complex_variable(dl[0]));
5750 GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
5751 "Data for coefficient should be a scalar");
5754 gmm::scaled(ccL, gmm::abs((*COEFF)[0])), vecl[0]);
5756 gmm::scaled(cB, gmm::abs((*COEFF)[0])), matl[0]);
5764 virtual std::string declare_volume_assembly_string
5765 (
const model &,
size_type,
const model::varnamelist &,
5766 const model::varnamelist &)
const {
5767 return std::string();
5770 constraint_brick(
bool penalized) {
5771 set_flags(penalized ?
"Constraint with penalization brick"
5772 :
"Constraint with multipliers brick",
5781 model_real_sparse_matrix &set_private_data_brick_real_matrix
5783 pbrick pbr = md.brick_pointer(indbrick);
5784 md.touch_brick(indbrick);
5785 have_private_data_brick *p =
dynamic_cast<have_private_data_brick *
>
5786 (
const_cast<virtual_brick *
>(pbr.get()));
5787 GMM_ASSERT1(p,
"Wrong type of brick");
5791 model_real_plain_vector &set_private_data_brick_real_rhs
5793 pbrick pbr = md.brick_pointer(indbrick);
5794 md.touch_brick(indbrick);
5795 have_private_data_brick *p =
dynamic_cast<have_private_data_brick *
>
5796 (
const_cast<virtual_brick *
>(pbr.get()));
5797 GMM_ASSERT1(p,
"Wrong type of brick");
5798 if (p->nameL.compare(
"") != 0) GMM_WARNING1(
"Rhs already set by data name");
5802 model_complex_sparse_matrix &set_private_data_brick_complex_matrix
5804 pbrick pbr = md.brick_pointer(indbrick);
5805 md.touch_brick(indbrick);
5806 have_private_data_brick *p =
dynamic_cast<have_private_data_brick *
>
5807 (
const_cast<virtual_brick *
>(pbr.get()));
5808 GMM_ASSERT1(p,
"Wrong type of brick");
5812 model_complex_plain_vector &set_private_data_brick_complex_rhs
5814 pbrick pbr = md.brick_pointer(indbrick);
5815 md.touch_brick(indbrick);
5816 have_private_data_brick *p =
dynamic_cast<have_private_data_brick *
>
5817 (
const_cast<virtual_brick *
>(pbr.get()));
5818 GMM_ASSERT1(p,
"Wrong type of brick");
5819 if (p->nameL.compare(
"") != 0) GMM_WARNING1(
"Rhs already set by data name");
5823 void set_private_data_rhs
5824 (model &md,
size_type indbrick,
const std::string &varname) {
5825 pbrick pbr = md.brick_pointer(indbrick);
5826 md.touch_brick(indbrick);
5827 have_private_data_brick *p =
dynamic_cast<have_private_data_brick *
>
5828 (
const_cast<virtual_brick *
>(pbr.get()));
5829 GMM_ASSERT1(p,
"Wrong type of brick");
5830 if (p->nameL.compare(varname) != 0) {
5831 model::varnamelist dl = md.datanamelist_of_brick(indbrick);
5832 if (p->nameL.compare(
"") == 0) dl.push_back(varname);
5833 else dl.back() = varname;
5834 md.change_data_of_brick(indbrick, dl);
5839 size_type add_constraint_with_penalization
5840 (model &md,
const std::string &varname, scalar_type penalisation_coeff) {
5841 std::string coeffname = md.new_name(
"penalization_on_" + varname);
5842 md.add_fixed_size_data(coeffname, 1);
5843 if (md.is_complex())
5844 md.set_complex_variable(coeffname)[0] = penalisation_coeff;
5846 md.set_real_variable(coeffname)[0] = penalisation_coeff;
5847 pbrick pbr = std::make_shared<constraint_brick>(
true);
5849 tl.push_back(model::term_description(varname, varname,
true));
5850 model::varnamelist vl(1, varname);
5851 model::varnamelist dl(1, coeffname);
5852 return md.add_brick(pbr, vl, dl, tl, model::mimlist(),
size_type(-1));
5855 size_type add_constraint_with_multipliers
5856 (model &md,
const std::string &varname,
const std::string &multname) {
5857 pbrick pbr = std::make_shared<constraint_brick>(
false);
5859 tl.push_back(model::term_description(multname, varname,
true));
5860 model::varnamelist vl(1, varname);
5861 vl.push_back(multname);
5862 model::varnamelist dl;
5863 return md.add_brick(pbr, vl, dl, tl, model::mimlist(),
size_type(-1));
5873 struct explicit_matrix_brick :
public have_private_data_brick {
5875 virtual void real_pre_assembly_in_serial(
const model &,
size_type,
5876 const model::varnamelist &vl,
5877 const model::varnamelist &dl,
5878 const model::mimlist &mims,
5879 model::real_matlist &matl,
5880 model::real_veclist &vecl,
5881 model::real_veclist &,
5883 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5884 "Explicit matrix has one and only one term");
5885 GMM_ASSERT1(mims.size() == 0,
"Explicit matrix need no mesh_im");
5886 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() == 0,
5887 "Wrong number of variables for explicit matrix brick");
5888 GMM_ASSERT1(gmm::mat_ncols(rB) == gmm::mat_ncols(matl[0]) &&
5889 gmm::mat_nrows(rB) == gmm::mat_nrows(matl[0]),
5890 "Explicit matrix brick dimension mismatch ("<<
5891 gmm::mat_ncols(rB)<<
"x"<<gmm::mat_nrows(rB)<<
") != ("<<
5892 gmm::mat_ncols(matl[0])<<
"x"<<gmm::mat_nrows(matl[0])<<
")");
5896 virtual void complex_pre_assembly_in_serial(
const model &,
size_type,
5897 const model::varnamelist &vl,
5898 const model::varnamelist &dl,
5899 const model::mimlist &mims,
5900 model::complex_matlist &matl,
5901 model::complex_veclist &vecl,
5902 model::complex_veclist &,
5904 build_version)
const {
5905 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5906 "Explicit matrix has one and only one term");
5907 GMM_ASSERT1(mims.size() == 0,
"Explicit matrix need no mesh_im");
5908 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() == 0,
5909 "Wrong number of variables for explicit matrix brick");
5913 virtual std::string declare_volume_assembly_string
5914 (
const model &,
size_type,
const model::varnamelist &,
5915 const model::varnamelist &)
const {
5916 return std::string();
5919 explicit_matrix_brick(
bool symmetric_,
bool coercive_) {
5920 set_flags(
"Explicit matrix brick",
5922 symmetric_ , coercive_ ,
5929 (model &md,
const std::string &varname1,
const std::string &varname2,
5930 bool issymmetric,
bool iscoercive) {
5931 pbrick pbr = std::make_shared<explicit_matrix_brick>(issymmetric,
5934 tl.push_back(model::term_description(varname1, varname2, issymmetric));
5935 model::varnamelist vl(1, varname1);
5936 vl.push_back(varname2);
5937 model::varnamelist dl;
5938 return md.add_brick(pbr, vl, dl, tl, model::mimlist(),
size_type(-1));
5947 struct explicit_rhs_brick :
public have_private_data_brick {
5949 virtual void real_pre_assembly_in_serial(
const model &,
size_type,
5950 const model::varnamelist &vl,
5951 const model::varnamelist &dl,
5952 const model::mimlist &mims,
5953 model::real_matlist &matl,
5954 model::real_veclist &vecl,
5955 model::real_veclist &,
5957 if (MPI_IS_MASTER()) {
5958 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5959 "Explicit rhs has one and only one term");
5960 GMM_ASSERT1(mims.size() == 0,
"Explicit rhs need no mesh_im");
5961 GMM_ASSERT1(vl.size() == 1 && dl.size() == 0,
5962 "Wrong number of variables for explicit rhs brick");
5967 virtual void complex_pre_assembly_in_serial(
const model &,
size_type,
5968 const model::varnamelist &vl,
5969 const model::varnamelist &dl,
5970 const model::mimlist &mims,
5971 model::complex_matlist &matl,
5972 model::complex_veclist &vecl,
5973 model::complex_veclist &,
5975 build_version)
const {
5976 if (MPI_IS_MASTER()) {
5977 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5978 "Explicit rhs has one and only one term");
5979 GMM_ASSERT1(mims.size() == 0,
"Explicit rhs need no mesh_im");
5980 GMM_ASSERT1(vl.size() == 1 && dl.size() == 0,
5981 "Wrong number of variables for explicit rhs brick");
5987 virtual std::string declare_volume_assembly_string
5988 (
const model &,
size_type,
const model::varnamelist &,
5989 const model::varnamelist &)
const {
5990 return std::string();
5993 explicit_rhs_brick() {
5994 set_flags(
"Explicit rhs brick",
6004 (model &md,
const std::string &varname) {
6005 pbrick pbr = std::make_shared<explicit_rhs_brick>();
6007 tl.push_back(model::term_description(varname));
6008 model::varnamelist vl(1, varname);
6009 model::varnamelist dl;
6010 return md.add_brick(pbr, vl, dl, tl, model::mimlist(),
size_type(-1));
6020 struct iso_lin_elasticity_new_brick :
public virtual_brick {
6022 std::string expr, dataname3;
6024 void asm_real_tangent_terms(
const model &md,
size_type ib,
6025 const model::varnamelist &vl,
6026 const model::varnamelist &dl,
6027 const model::mimlist &mims,
6028 model::real_matlist &matl,
6029 model::real_veclist &vecl,
6030 model::real_veclist &,
6032 build_version version)
const override {
6033 GMM_ASSERT1(vl.size() == 1,
"Linearized isotropic elasticity brick "
6034 "has one and only one variable");
6035 GMM_ASSERT1(matl.size() == 1,
"Linearized isotropic elasticity brick "
6036 "has one and only one term");
6037 GMM_ASSERT1(mims.size() == 1,
"Linearized isotropic elasticity brick "
6038 "needs one and only one mesh_im");
6040 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0);
6041 for (
size_type i = 0; i < dl.size(); ++i) {
6042 recompute_matrix = recompute_matrix ||
6043 md.is_var_newer_than_brick(dl[i], ib);
6046 if (recompute_matrix) {
6048 ga_workspace workspace(md, ga_workspace::inherit::ALL);
6049 workspace.add_expression(expr, *(mims[0]), region);
6050 GMM_TRACE2(name <<
": generic matrix assembly");
6051 workspace.assembly(2);
6052 scalar_type
alpha = scalar_type(1)
6053 / (workspace.factor_of_variable(vl[0]));
6054 const auto &R=workspace.assembled_matrix();
6055 gmm::sub_interval I = workspace.interval_of_variable(vl[0]);
6056 gmm::copy(gmm::scaled(gmm::sub_matrix(R, I, I), alpha),
6060 if (dataname3.size()) {
6066 gmm::scaled(md.real_variable(dataname3), scalar_type(-1)),
6072 void real_post_assembly_in_serial(
const model &md,
size_type ib,
6073 const model::varnamelist &,
6074 const model::varnamelist &,
6075 const model::mimlist &,
6076 model::real_matlist &,
6077 model::real_veclist &vecl,
6078 model::real_veclist &,
6080 build_version)
const override {
6081 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
6085 virtual std::string declare_volume_assembly_string
6086 (
const model &,
size_type,
const model::varnamelist &,
6087 const model::varnamelist &)
const {
6091 iso_lin_elasticity_new_brick(
const std::string &expr_,
6092 const std::string &dataname3_) {
6093 expr = expr_; dataname3 = dataname3_;
6094 set_flags(
"Linearized isotropic elasticity",
true ,
6103 (
model &md,
const mesh_im &mim,
const std::string &varname,
6104 const std::string &dataexpr1,
const std::string &dataexpr2,
6105 size_type region,
const std::string &dataname3) {
6106 std::string test_varname
6107 =
"Test_" + sup_previous_and_dot_to_varname(varname);
6109 std::string expr1 =
"((("+dataexpr1+
")*(Div_"+varname+
"-Div_"+dataname3
6110 +
"))*Id(meshdim)+(2*("+dataexpr2+
"))*(Sym(Grad_"+varname
6111 +
")-Sym(Grad_"+dataname3+
"))):Grad_" +test_varname;
6112 std::string expr2 =
"(Div_"+varname+
"*(("+dataexpr1+
")*Id(meshdim))"
6113 +
"+(2*("+dataexpr2+
"))*Sym(Grad_"+varname+
")):Grad_"+test_varname;
6116 model::varnamelist vl, dl;
6118 ga_workspace workspace(md, ga_workspace::inherit::ALL);
6119 workspace.add_expression(expr2, mim, region);
6120 model::varnamelist vl_test1, vl_test2;
6121 is_lin = workspace.used_variables(vl, vl_test1, vl_test2, dl, 2);
6124 pbrick pbr = std::make_shared<iso_lin_elasticity_new_brick>
6127 tl.push_back(model::term_description(varname,
6128 sup_previous_and_dot_to_varname(varname),
true));
6129 if (dataname3.size()) dl.push_back(dataname3);
6130 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
6133 (md, mim, dataname3.size() ? expr1 : expr2, region,
false,
false,
6134 "Linearized isotropic elasticity (with nonlinear dependance)");
6139 (
model &md,
const mesh_im &mim,
const std::string &varname,
6140 const std::string &data_E,
const std::string &data_nu,
6142 std::string test_varname
6143 =
"Test_" + sup_previous_and_dot_to_varname(varname);
6145 std::string mu =
"(("+data_E+
")/(2*(1+("+data_nu+
"))))";
6146 std::string lambda =
"(("+data_E+
")*("+data_nu+
")/((1+("+data_nu+
"))*(1-2*("
6148 std::string expr = lambda+
"*Div_"+varname+
"*Div_"+test_varname
6149 +
"+"+mu+
"*(Grad_"+varname+
"+Grad_"+varname+
"'):Grad_"+test_varname;
6153 ga_workspace workspace(md, ga_workspace::inherit::ALL);
6154 workspace.add_expression(expr, mim, region);
6155 is_lin = workspace.is_linear(2);
6159 "Linearized isotropic elasticity");
6162 (md, mim, expr, region,
false,
false,
6163 "Linearized isotropic elasticity (with nonlinear dependance)");
6168 (
model &md,
const mesh_im &mim,
const std::string &varname,
6169 const std::string &data_E,
const std::string &data_nu,
6171 std::string test_varname
6172 =
"Test_" + sup_previous_and_dot_to_varname(varname);
6175 GMM_ASSERT1(mfu,
"The variable should be a fem variable");
6178 std::string mu =
"(("+data_E+
")/(2*(1+("+data_nu+
"))))";
6179 std::string lambda =
"(("+data_E+
")*("+data_nu+
")/((1+("+data_nu+
"))*(1-2*("
6182 lambda =
"(("+data_E+
")*("+data_nu+
")/((1-sqr("+data_nu+
"))))";
6183 std::string expr = lambda+
"*Div_"+varname+
"*Div_"+test_varname
6184 +
"+"+mu+
"*(Grad_"+varname+
"+Grad_"+varname+
"'):Grad_"+test_varname;
6188 ga_workspace workspace(md, ga_workspace::inherit::ALL);
6189 workspace.add_expression(expr, mim, region);
6190 is_lin = workspace.is_linear(2);
6194 "Linearized isotropic elasticity");
6197 (md, mim, expr, region,
false,
false,
6198 "Linearized isotropic elasticity (with nonlinear dependance)");
6203 void compute_isotropic_linearized_Von_Mises_or_Tresca
6204 (model &md,
const std::string &varname,
const std::string &data_lambda,
6205 const std::string &data_mu,
const mesh_fem &mf_vm,
6206 model_real_plain_vector &VM,
bool tresca) {
6209 const mesh_fem &mf_u = md.mesh_fem_of_variable(varname);
6210 const mesh_fem *mf_lambda = md.pmesh_fem_of_variable(data_lambda);
6211 const model_real_plain_vector *lambda=&(md.real_variable(data_lambda));
6212 const mesh_fem *mf_mu = md.pmesh_fem_of_variable(data_mu);
6213 const model_real_plain_vector *mu = &(md.real_variable(data_mu));
6216 if (mf_lambda) sl = sl * mf_lambda->get_qdim() / mf_lambda->nb_dof();
6218 if (mf_mu) sm = sm * mf_mu->get_qdim() / mf_mu->nb_dof();
6220 GMM_ASSERT1(sl == 1 && sm == 1,
"Bad format for Lame coefficients");
6221 GMM_ASSERT1(mf_lambda == mf_mu,
6222 "The two Lame coefficients should be described on the same "
6223 "finite element method.");
6227 md.real_variable(varname), VM,
6228 *mf_lambda, *lambda,
6233 model_real_plain_vector LAMBDA(mf_lambda->nb_dof(), (*lambda)[0]);
6234 model_real_plain_vector MU(mf_lambda->nb_dof(), (*mu)[0]);
6236 md.real_variable(varname), VM,
6245 std::string sigma_d=
"("+data_mu+
")*(Grad_"+varname+
"+Grad_"+varname+
"')";
6246 std::string expr =
"sqrt(3/2)*Norm(Deviator("+sigma_d+
"))";
6247 ga_interpolation_Lagrange_fem(md, expr, mf_vm, VM);
6253 (
model &md,
const std::string &varname,
const std::string &data_E,
6254 const std::string &data_nu,
const mesh_fem &mf_vm,
6255 model_real_plain_vector &VM) {
6259 std::string mu =
"(("+data_E+
")/(2*(1+("+data_nu+
"))))";
6262 std::string sigma_d = mu+
"*(Grad_"+varname+
"+Grad_"+varname+
"')";
6263 std::string expr =
"sqrt(3/2)*Norm(Deviator("+sigma_d+
"))";
6264 ga_interpolation_Lagrange_fem(md, expr, mf_vm, VM);
6268 (
model &md,
const std::string &varname,
const std::string &data_E,
6269 const std::string &data_nu,
const mesh_fem &mf_vm,
6270 model_real_plain_vector &VM) {
6279 std::string mu =
"(("+data_E+
")/(2*(1+("+data_nu+
"))))";
6282 std::string sigma_d = mu+
"*(Grad_"+varname+
"+Grad_"+varname+
"')";
6283 std::string expr =
"sqrt(3/2)*Norm(Deviator("+sigma_d+
"))";
6284 ga_interpolation_Lagrange_fem(md, expr, mf_vm, VM);
6294 struct linear_incompressibility_brick :
public virtual_brick {
6296 virtual void asm_real_tangent_terms(
const model &md,
size_type ,
6297 const model::varnamelist &vl,
6298 const model::varnamelist &dl,
6299 const model::mimlist &mims,
6300 model::real_matlist &matl,
6301 model::real_veclist &,
6302 model::real_veclist &,
6304 build_version)
const {
6306 GMM_ASSERT1((matl.size() == 1 && dl.size() == 0)
6307 || (matl.size() == 2 && dl.size() == 1),
6308 "Wrong term and/or data number for Linear incompressibility "
6310 GMM_ASSERT1(mims.size() == 1,
"Linear incompressibility brick need one "
6311 "and only one mesh_im");
6312 GMM_ASSERT1(vl.size() == 2,
"Wrong number of variables for linear "
6313 "incompressibility brick");
6315 bool penalized = (dl.size() == 1);
6316 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6317 const mesh_fem &mf_p = md.mesh_fem_of_variable(vl[1]);
6318 const mesh_im &mim = *mims[0];
6319 const model_real_plain_vector *COEFF = 0;
6320 const mesh_fem *mf_data = 0;
6323 COEFF = &(md.real_variable(dl[0]));
6324 mf_data = md.pmesh_fem_of_variable(dl[0]);
6326 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
6327 GMM_ASSERT1(s == 1,
"Bad format for the penalization parameter");
6330 mesh_region rg(region);
6331 mim.linked_mesh().intersect_with_mpi_region(rg);
6333 GMM_TRACE2(
"Stokes term assembly");
6341 gmm::scale(matl[1], scalar_type(-1));
6345 gmm::scale(matl[1], -(*COEFF)[0]);
6352 virtual void real_post_assembly_in_serial(
const model &,
size_type,
6353 const model::varnamelist &,
6354 const model::varnamelist &,
6355 const model::mimlist &,
6356 model::real_matlist &,
6357 model::real_veclist &,
6358 model::real_veclist &,
6360 build_version)
const
6364 linear_incompressibility_brick() {
6365 set_flags(
"Linear incompressibility brick",
6374 (
model &md,
const mesh_im &mim,
const std::string &varname,
6375 const std::string &multname,
size_type region,
6376 const std::string &dataexpr) {
6378 pbrick pbr = std::make_shared<linear_incompressibility_brick>();
6380 tl.push_back(model::term_description(multname, varname,
true));
6381 model::varnamelist vl(1, varname);
6382 vl.push_back(multname);
6383 model::varnamelist dl;
6384 if (dataname.size()) {
6385 dl.push_back(dataexpr);
6386 tl.push_back(model::term_description(multname, multname,
true));
6388 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
6390 std::string test_varname
6391 =
"Test_" + sup_previous_and_dot_to_varname(varname);
6392 std::string test_multname
6393 =
"Test_" + sup_previous_and_dot_to_varname(multname);
6395 if (dataexpr.size())
6396 expr =
"-"+multname+
"*Div_"+test_varname +
"-"+test_multname
6397 +
"*Div_"+varname+
"+(("+dataexpr+
")*"+multname+
")*"+test_multname;
6399 expr =
"-"+multname+
"*Div_"+test_varname +
"-"+test_multname
6402 "Linear incompressibility",
true);
6405 (md, mim, expr, region,
false,
false,
6406 "Linear incompressibility (with nonlinear dependance)");
6419 struct mass_brick :
public virtual_brick {
6421 virtual void asm_real_tangent_terms(
const model &md,
size_type,
6422 const model::varnamelist &vl,
6423 const model::varnamelist &dl,
6424 const model::mimlist &mims,
6425 model::real_matlist &matl,
6426 model::real_veclist &,
6427 model::real_veclist &,
6429 build_version)
const {
6430 GMM_ASSERT1(matl.size() == 1,
6431 "Mass brick has one and only one term");
6432 GMM_ASSERT1(mims.size() == 1,
6433 "Mass brick need one and only one mesh_im");
6434 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
6435 "Wrong number of variables for mass brick");
6437 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6438 const mesh &m = mf_u.linked_mesh();
6439 const mesh_im &mim = *mims[0];
6440 mesh_region rg(region);
6441 m.intersect_with_mpi_region(rg);
6443 const mesh_fem *mf_rho = 0;
6444 const model_real_plain_vector *rho = 0;
6447 mf_rho = md.pmesh_fem_of_variable(dl[0]);
6448 rho = &(md.real_variable(dl[0]));
6450 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6451 GMM_ASSERT1(sl == 1,
"Bad format of mass brick coefficient");
6454 GMM_TRACE2(
"Mass matrix assembly");
6456 if (dl.size() && mf_rho) {
6460 if (dl.size()) gmm::scale(matl[0], (*rho)[0]);
6464 virtual void asm_complex_tangent_terms(
const model &md,
size_type,
6465 const model::varnamelist &vl,
6466 const model::varnamelist &dl,
6467 const model::mimlist &mims,
6468 model::complex_matlist &matl,
6469 model::complex_veclist &,
6470 model::complex_veclist &,
6472 build_version)
const {
6473 GMM_ASSERT1(matl.size() == 1,
6474 "Mass brick has one and only one term");
6475 GMM_ASSERT1(mims.size() == 1,
6476 "Mass brick need one and only one mesh_im");
6477 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
6478 "Wrong number of variables for mass brick");
6480 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6481 const mesh &m = mf_u.linked_mesh();
6482 const mesh_im &mim = *mims[0];
6483 mesh_region rg(region);
6484 m.intersect_with_mpi_region(rg);
6486 const mesh_fem *mf_rho = 0;
6487 const model_complex_plain_vector *rho = 0;
6490 mf_rho = md.pmesh_fem_of_variable(dl[0]);
6491 rho = &(md.complex_variable(dl[0]));
6493 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6494 GMM_ASSERT1(sl == 1,
"Bad format of mass brick coefficient");
6497 GMM_TRACE2(
"Mass matrix assembly");
6499 if (dl.size() && mf_rho) {
6503 if (dl.size()) gmm::scale(matl[0], (*rho)[0]);
6507 virtual std::string declare_volume_assembly_string
6508 (
const model &,
size_type,
const model::varnamelist &,
6509 const model::varnamelist &)
const {
6510 return std::string();
6514 set_flags(
"Mass brick",
true ,
6523 (
model &md,
const mesh_im &mim,
const std::string &varname,
6524 const std::string &dataexpr_rho,
size_type region) {
6526 pbrick pbr = std::make_shared<mass_brick>();
6528 tl.push_back(model::term_description(varname, varname,
true));
6529 model::varnamelist dl;
6530 if (dataexpr_rho.size())
6531 dl.push_back(dataexpr_rho);
6532 return md.
add_brick(pbr, model::varnamelist(1, varname), dl, tl,
6533 model::mimlist(1, &mim), region);
6535 std::string test_varname
6536 =
"Test_" + sup_previous_and_dot_to_varname(varname);
6538 if (dataexpr_rho.size())
6539 expr =
"(("+dataexpr_rho+
")*"+varname+
")."+test_varname;
6541 expr = varname+
"."+test_varname;
6543 "Mass matrix",
true);
6546 "Mass matrix (nonlinear)");
6557 struct lumped_mass_for_first_order_brick :
public virtual_brick {
6559 virtual void asm_real_tangent_terms(
const model &md,
size_type,
6560 const model::varnamelist &vl,
6561 const model::varnamelist &dl,
6562 const model::mimlist &mims,
6563 model::real_matlist &matl,
6564 model::real_veclist &,
6565 model::real_veclist &,
6567 build_version)
const {
6568 GMM_ASSERT1(matl.size() == 1,
6569 "Lumped Mass brick has one and only one term");
6570 GMM_ASSERT1(mims.size() == 1,
6571 "Lumped Mass brick needs one and only one mesh_im");
6572 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
6573 "Wrong number of variables for lumped mass brick");
6575 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6576 const mesh &m = mf_u.linked_mesh();
6577 const mesh_im &mim = *mims[0];
6578 mesh_region rg(region);
6579 m.intersect_with_mpi_region(rg);
6581 const mesh_fem *mf_rho = 0;
6582 const model_real_plain_vector *rho = 0;
6585 mf_rho = md.pmesh_fem_of_variable(dl[0]);
6586 rho = &(md.real_variable(dl[0]));
6588 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6589 GMM_ASSERT1(sl == 1,
"Bad format of mass brick coefficient");
6592 GMM_TRACE2(
"Lumped mass matrix assembly (please check that integration is 1st order.)");
6594 if (dl.size() && mf_rho) {
6598 if (dl.size()) gmm::scale(matl[0], (*rho)[0]);
6603 lumped_mass_for_first_order_brick() {
6604 set_flags(
"Lumped mass brick",
true ,
6613 (
model & md,
const mesh_im &mim,
const std::string &varname,
6614 const std::string &dataexpr_rho,
size_type region) {
6615 pbrick pbr = std::make_shared<lumped_mass_for_first_order_brick>();
6617 tl.push_back(model::term_description(varname, varname,
true));
6618 model::varnamelist dl;
6619 if (dataexpr_rho.size())
6620 dl.push_back(dataexpr_rho);
6621 return md.
add_brick(pbr, model::varnamelist(1, varname), dl, tl,
6622 model::mimlist(1, &mim), region);
6638 struct basic_d_on_dt_brick :
public virtual_brick {
6640 virtual void asm_real_tangent_terms(
const model &md,
size_type ib,
6641 const model::varnamelist &vl,
6642 const model::varnamelist &dl,
6643 const model::mimlist &mims,
6644 model::real_matlist &matl,
6645 model::real_veclist &vecl,
6646 model::real_veclist &,
6648 build_version version)
const {
6649 GMM_ASSERT1(matl.size() == 1,
6650 "Basic d/dt brick has one and only one term");
6651 GMM_ASSERT1(mims.size() == 1,
6652 "Basic d/dt brick need one and only one mesh_im");
6653 GMM_ASSERT1(vl.size() == 1 && dl.size() >= 2 && dl.size() <= 3,
6654 "Wrong number of variables for basic d/dt brick");
6658 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
6659 || (md.is_var_newer_than_brick(dl[1], ib));
6661 recompute_matrix = recompute_matrix ||
6662 md.is_var_newer_than_brick(dl[2], ib);
6665 if (recompute_matrix) {
6666 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6667 const mesh &m = mf_u.linked_mesh();
6668 const mesh_im &mim = *mims[0];
6669 mesh_region rg(region);
6670 m.intersect_with_mpi_region(rg);
6672 const model_real_plain_vector &dt = md.real_variable(dl[1]);
6673 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for time step");
6675 const mesh_fem *mf_rho = 0;
6676 const model_real_plain_vector *rho = 0;
6678 if (dl.size() > 2) {
6679 mf_rho = md.pmesh_fem_of_variable(dl[2]);
6680 rho = &(md.real_variable(dl[2]));
6682 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6683 GMM_ASSERT1(sl == 1,
"Bad format for density");
6686 GMM_TRACE2(
"Mass matrix assembly for d_on_dt brick");
6687 if (dl.size() > 2 && mf_rho) {
6690 gmm::scale(matl[0], scalar_type(1) / dt[0]);
6694 if (dl.size() > 2) gmm::scale(matl[0], (*rho)[0] / dt[0]);
6695 else gmm::scale(matl[0], scalar_type(1) / dt[0]);
6698 gmm::mult(matl[0], md.real_variable(dl[0], 1), vecl[0]);
6701 virtual void asm_complex_tangent_terms(
const model &md,
size_type ib,
6702 const model::varnamelist &vl,
6703 const model::varnamelist &dl,
6704 const model::mimlist &mims,
6705 model::complex_matlist &matl,
6706 model::complex_veclist &vecl,
6707 model::complex_veclist &,
6709 build_version version)
const {
6710 GMM_ASSERT1(matl.size() == 1,
6711 "Basic d/dt brick has one and only one term");
6712 GMM_ASSERT1(mims.size() == 1,
6713 "Basic d/dt brick need one and only one mesh_im");
6714 GMM_ASSERT1(vl.size() == 1 && dl.size() >= 2 && dl.size() <= 3,
6715 "Wrong number of variables for basic d/dt brick");
6719 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
6720 || (md.is_var_newer_than_brick(dl[1], ib));
6722 recompute_matrix = recompute_matrix ||
6723 md.is_var_newer_than_brick(dl[2], ib);
6725 if (recompute_matrix) {
6726 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6727 const mesh &m = mf_u.linked_mesh();
6728 const mesh_im &mim = *mims[0];
6730 mesh_region rg(region);
6731 m.intersect_with_mpi_region(rg);
6733 const model_complex_plain_vector &dt = md.complex_variable(dl[1]);
6734 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for time step");
6736 const mesh_fem *mf_rho = 0;
6737 const model_complex_plain_vector *rho = 0;
6739 if (dl.size() > 2) {
6740 mf_rho = md.pmesh_fem_of_variable(dl[2]);
6741 rho = &(md.complex_variable(dl[2]));
6743 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6744 GMM_ASSERT1(sl == 1,
"Bad format for density");
6747 GMM_TRACE2(
"Mass matrix assembly for d_on_dt brick");
6748 if (dl.size() > 2 && mf_rho) {
6751 gmm::scale(matl[0], scalar_type(1) / dt[0]);
6755 if (dl.size() > 2) gmm::scale(matl[0], (*rho)[0] / dt[0]);
6756 else gmm::scale(matl[0], scalar_type(1) / dt[0]);
6759 gmm::mult(matl[0], md.complex_variable(dl[0], 1), vecl[0]);
6762 virtual std::string declare_volume_assembly_string
6763 (
const model &,
size_type,
const model::varnamelist &,
6764 const model::varnamelist &)
const {
6765 return std::string();
6768 basic_d_on_dt_brick() {
6769 set_flags(
"Basic d/dt brick",
true ,
6778 (
model &md,
const mesh_im &mim,
const std::string &varname,
6779 const std::string &dataname_dt,
const std::string &dataname_rho,
6781 pbrick pbr = std::make_shared<basic_d_on_dt_brick>();
6783 tl.push_back(model::term_description(varname, varname,
true));
6784 model::varnamelist dl(1, varname);
6785 dl.push_back(dataname_dt);
6786 if (dataname_rho.size())
6787 dl.push_back(dataname_rho);
6788 return md.
add_brick(pbr, model::varnamelist(1, varname), dl, tl,
6789 model::mimlist(1, &mim), region);
6800 struct basic_d2_on_dt2_brick :
public virtual_brick {
6802 mutable scalar_type old_alphadt2;
6804 virtual void asm_real_tangent_terms(
const model &md,
size_type ib,
6805 const model::varnamelist &vl,
6806 const model::varnamelist &dl,
6807 const model::mimlist &mims,
6808 model::real_matlist &matl,
6809 model::real_veclist &vecl,
6810 model::real_veclist &,
6812 build_version version)
const {
6813 GMM_ASSERT1(matl.size() == 1,
6814 "Basic d2/dt2 brick has one and only one term");
6815 GMM_ASSERT1(mims.size() == 1,
6816 "Basic d2/dt2 brick need one and only one mesh_im");
6817 GMM_ASSERT1(vl.size() == 1 && dl.size() >= 4 && dl.size() <= 5,
6818 "Wrong number of variables for basic d2/dt2 brick");
6820 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0);
6822 recompute_matrix = recompute_matrix || md.is_var_newer_than_brick(dl[2], ib);
6823 if (dl.size() > 4) recompute_matrix || md.is_var_newer_than_brick(dl[4], ib);
6825 const model_real_plain_vector &dt = md.real_variable(dl[2]);
6826 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for time step");
6827 const model_real_plain_vector &alpha = md.real_variable(dl[3]);
6828 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for parameter alpha");
6829 scalar_type alphadt2 = gmm::sqr(dt[0]) * alpha[0];
6831 if (!recompute_matrix && alphadt2 != old_alphadt2)
6832 gmm::scale(matl[0], old_alphadt2/alphadt2);
6833 old_alphadt2 = alphadt2;
6835 if (recompute_matrix) {
6836 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6837 const mesh &m = mf_u.linked_mesh();
6838 const mesh_im &mim = *mims[0];
6839 mesh_region rg(region);
6840 m.intersect_with_mpi_region(rg);
6842 const mesh_fem *mf_rho = 0;
6843 const model_real_plain_vector *rho = 0;
6845 if (dl.size() > 4) {
6846 mf_rho = md.pmesh_fem_of_variable(dl[4]);
6847 rho = &(md.real_variable(dl[4]));
6849 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6850 GMM_ASSERT1(sl == 1,
"Bad format for density");
6853 GMM_TRACE2(
"Mass matrix assembly for d2_on_dt2 brick");
6854 if (dl.size() > 4 && mf_rho) {
6857 gmm::scale(matl[0], scalar_type(1) / alphadt2);
6861 if (dl.size() > 4) gmm::scale(matl[0], (*rho)[0] / alphadt2);
6862 else gmm::scale(matl[0], scalar_type(1) / alphadt2);
6865 gmm::mult(matl[0], md.real_variable(dl[0], 1), vecl[0]);
6866 gmm::mult_add(matl[0], gmm::scaled(md.real_variable(dl[1], 1), dt[0]),
6870 virtual void asm_complex_tangent_terms(
const model &md,
size_type ib,
6871 const model::varnamelist &vl,
6872 const model::varnamelist &dl,
6873 const model::mimlist &mims,
6874 model::complex_matlist &matl,
6875 model::complex_veclist &vecl,
6876 model::complex_veclist &,
6878 build_version version)
const {
6879 GMM_ASSERT1(matl.size() == 1,
6880 "Basic d2/dt2 brick has one and only one term");
6881 GMM_ASSERT1(mims.size() == 1,
6882 "Basic d2/dt2 brick need one and only one mesh_im");
6883 GMM_ASSERT1(vl.size() == 1 && dl.size() >= 4 && dl.size() <= 5,
6884 "Wrong number of variables for basic d2/dt2 brick");
6886 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0);
6888 recompute_matrix = recompute_matrix || md.is_var_newer_than_brick(dl[2], ib);
6889 if (dl.size() > 4) recompute_matrix || md.is_var_newer_than_brick(dl[4], ib);
6892 const model_complex_plain_vector &dt = md.complex_variable(dl[2]);
6893 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for time step");
6894 const model_complex_plain_vector &
alpha = md.complex_variable(dl[3]);
6895 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for parameter alpha");
6896 scalar_type alphadt2 = gmm::real(gmm::sqr(dt[0]) * alpha[0]);
6898 if (!recompute_matrix && alphadt2 != old_alphadt2)
6899 gmm::scale(matl[0], old_alphadt2/alphadt2);
6900 old_alphadt2 = alphadt2;
6902 if (recompute_matrix) {
6903 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6904 const mesh &m = mf_u.linked_mesh();
6905 const mesh_im &mim = *mims[0];
6906 mesh_region rg(region);
6907 m.intersect_with_mpi_region(rg);
6909 const mesh_fem *mf_rho = 0;
6910 const model_complex_plain_vector *rho = 0;
6912 if (dl.size() > 4) {
6913 mf_rho = md.pmesh_fem_of_variable(dl[4]);
6914 rho = &(md.complex_variable(dl[4]));
6916 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6917 GMM_ASSERT1(sl == 1,
"Bad format for density");
6920 GMM_TRACE2(
"Mass matrix assembly for d2_on_dt2 brick");
6921 if (dl.size() > 4 && mf_rho) {
6924 gmm::scale(matl[0], scalar_type(1) / alphadt2);
6928 if (dl.size() > 4) gmm::scale(matl[0], (*rho)[0] / alphadt2);
6929 else gmm::scale(matl[0], scalar_type(1) / alphadt2);
6932 gmm::mult(matl[0], md.complex_variable(dl[0], 1), vecl[0]);
6933 gmm::mult_add(matl[0], gmm::scaled(md.complex_variable(dl[1], 1), dt[0]),
6937 virtual std::string declare_volume_assembly_string
6938 (
const model &,
size_type,
const model::varnamelist &,
6939 const model::varnamelist &)
const {
6940 return std::string();
6943 basic_d2_on_dt2_brick() {
6944 set_flags(
"Basic d2/dt2 brick",
true ,
6953 (
model &md,
const mesh_im &mim,
const std::string &varnameU,
6954 const std::string &datanameV,
6955 const std::string &dataname_dt,
6956 const std::string &dataname_alpha,
6957 const std::string &dataname_rho,
6959 pbrick pbr = std::make_shared<basic_d2_on_dt2_brick>();
6961 tl.push_back(model::term_description(varnameU, varnameU,
true));
6962 model::varnamelist dl(1, varnameU);
6963 dl.push_back(datanameV);
6964 dl.push_back(dataname_dt);
6965 dl.push_back(dataname_alpha);
6966 if (dataname_rho.size())
6967 dl.push_back(dataname_rho);
6968 return md.
add_brick(pbr, model::varnamelist(1, varnameU), dl, tl,
6969 model::mimlist(1, &mim), region);
6988 void theta_method_dispatcher::set_dispatch_coeff(
const model &md,
size_type ib)
const {
6990 if (md.is_complex())
6991 theta = gmm::real(md.complex_variable(param_names[0])[0]);
6993 theta = md.real_variable(param_names[0])[0];
6995 md.matrix_coeff_of_brick(ib) = theta;
6997 md.rhs_coeffs_of_brick(ib)[0] = theta;
6999 md.rhs_coeffs_of_brick(ib)[1] = (scalar_type(1) - theta);
7003 void theta_method_dispatcher::next_real_iter
7004 (
const model &md,
size_type ib,
const model::varnamelist &vl,
7005 const model::varnamelist &dl, model::real_matlist &matl,
7006 std::vector<model::real_veclist> &vectl,
7007 std::vector<model::real_veclist> &vectl_sym,
bool first_iter)
const {
7008 next_iter(md, ib, vl, dl, matl, vectl, vectl_sym, first_iter);
7011 void theta_method_dispatcher::next_complex_iter
7012 (
const model &md,
size_type ib,
const model::varnamelist &vl,
7013 const model::varnamelist &dl,
7014 model::complex_matlist &matl,
7015 std::vector<model::complex_veclist> &vectl,
7016 std::vector<model::complex_veclist> &vectl_sym,
7017 bool first_iter)
const {
7018 next_iter(md, ib, vl, dl, matl, vectl, vectl_sym, first_iter);
7021 void theta_method_dispatcher::asm_real_tangent_terms
7022 (
const model &md,
size_type ib, model::real_matlist &,
7023 std::vector<model::real_veclist> &,
7024 std::vector<model::real_veclist> &,
7025 build_version version)
const
7026 { md.brick_call(ib, version, 0); }
7028 void theta_method_dispatcher::asm_complex_tangent_terms
7029 (
const model &md,
size_type ib, model::complex_matlist &,
7030 std::vector<model::complex_veclist> &,
7031 std::vector<model::complex_veclist> &,
7032 build_version version)
const
7033 { md.brick_call(ib, version, 0); }
7035 theta_method_dispatcher::theta_method_dispatcher(
const std::string &THETA)
7036 : virtual_dispatcher(2) {
7037 param_names.push_back(THETA);
7041 (
model &md, dal::bit_vector ibricks,
const std::string &THETA) {
7042 pdispatcher pdispatch = std::make_shared<theta_method_dispatcher>(THETA);
7043 for (dal::bv_visitor i(ibricks); !i.finished(); ++i)
7048 (
model &md,
const std::string &U,
const std::string &V,
7049 const std::string &pdt,
const std::string &ptheta) {
7055 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for time step");
7057 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for parameter theta");
7060 scalar_type(1) - scalar_type(1) / theta[0]),
7063 scalar_type(1) / (theta[0]*dt[0])),
7066 -scalar_type(1) / (theta[0]*dt[0])),
7070 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for time step");
7071 const model_real_plain_vector &theta = md.
real_variable(ptheta);
7072 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for parameter theta");
7075 scalar_type(1) - scalar_type(1) / theta[0]),
7078 scalar_type(1) / (theta[0]*dt[0])),
7081 -scalar_type(1) / (theta[0]*dt[0])),
7093 (
model &md,
size_type id2dt2b,
const std::string &U,
const std::string &V,
7094 const std::string &pdt,
const std::string &ptwobeta,
7095 const std::string &pgamma) {
7105 if (twobeta != gamma) {
7107 md.set_dispatch_coeff();
7111 md.
assembly(model::BUILD_RHS_WITH_LIN);
7114 model_complex_plain_vector W(nbdof), RHS(nbdof);
7115 gmm::copy(gmm::sub_vector(md.
complex_rhs(), md.interval_of_variable(U)),
7120 gmm::cg(md.linear_complex_matrix_term(id2dt2b, 0),
7121 W, RHS, gmm::identity_matrix(), iter);
7122 GMM_ASSERT1(iter.converged(),
"Velocity not well computed");
7124 gmm::scaled(W, complex_type(1)/(twobeta*dt)),
7128 if (twobeta != gamma) {
7130 md.set_dispatch_coeff();
7134 GMM_ASSERT1(
false,
"to be done");
7143 if (twobeta != gamma) {
7145 md.set_dispatch_coeff();
7149 md.
assembly(model::BUILD_RHS_WITH_LIN);
7152 model_real_plain_vector W(nbdof), RHS(nbdof);
7153 gmm::copy(gmm::sub_vector(md.
real_rhs(), md.interval_of_variable(U)),
7158 gmm::cg(md.linear_real_matrix_term(id2dt2b, 0),
7159 W, RHS, gmm::identity_matrix(), iter);
7160 GMM_ASSERT1(iter.converged(),
"Velocity not well computed");
7162 gmm::scaled(W, scalar_type(1)/(twobeta*dt)),
7166 if (twobeta != gamma) {
7168 md.set_dispatch_coeff();
7183 class midpoint_dispatcher :
public virtual_dispatcher {
7185 gmm::uint64_type id_num;
7189 typedef model::build_version build_version;
7191 void set_dispatch_coeff(
const model &md,
size_type ib)
const {
7192 md.matrix_coeff_of_brick(ib) = scalar_type(1)/scalar_type(2);
7193 md.rhs_coeffs_of_brick(ib)[0] = scalar_type(1);
7194 md.rhs_coeffs_of_brick(ib)[1] = scalar_type(1)/scalar_type(2);
7197 template <
typename MATLIST,
typename VECTLIST>
7198 inline void next_iter(
const model &md,
size_type ib,
7199 const model::varnamelist &vl,
7200 const model::varnamelist &dl,
7202 VECTLIST &vectl, VECTLIST &vectl_sym,
7203 bool first_iter)
const {
7205 pbrick pbr = md.brick_pointer(ib);
7209 if (!(pbr->is_linear()))
7210 md.add_temporaries(vl, id_num);
7211 md.add_temporaries(dl, id_num);
7216 if (pbr->is_linear()) {
7219 if (first_iter) md.update_brick(ib, model::BUILD_RHS);
7222 md.linear_brick_add_to_rhs(ib, 1, 0);
7227 (
const model &md,
size_type ib,
const model::varnamelist &vl,
7228 const model::varnamelist &dl, model::real_matlist &matl,
7229 std::vector<model::real_veclist> &vectl,
7230 std::vector<model::real_veclist> &vectl_sym,
bool first_iter)
const {
7231 next_iter(md, ib, vl, dl, matl, vectl, vectl_sym, first_iter);
7234 void next_complex_iter
7235 (
const model &md,
size_type ib,
const model::varnamelist &vl,
7236 const model::varnamelist &dl,
7237 model::complex_matlist &matl,
7238 std::vector<model::complex_veclist> &vectl,
7239 std::vector<model::complex_veclist> &vectl_sym,
7240 bool first_iter)
const {
7241 next_iter(md, ib, vl, dl, matl, vectl, vectl_sym, first_iter);
7244 void asm_real_tangent_terms
7245 (
const model &md,
size_type ib, model::real_matlist &,
7246 std::vector<model::real_veclist> &vectl,
7247 std::vector<model::real_veclist> &vectl_sym,
7248 build_version version)
const {
7250 scalar_type half = scalar_type(1)/scalar_type(2);
7251 pbrick pbr = md.brick_pointer(ib);
7254 const model::varnamelist &vl = md.varnamelist_of_brick(ib);
7255 const model::varnamelist &dl = md.datanamelist_of_brick(ib);
7257 if (!(pbr->is_linear())) {
7258 for (
size_type i = 0; i < vl.size(); ++i) {
7259 bool is_uptodate = md.temporary_uptodate(vl[i], id_num, ind);
7260 if (!is_uptodate && ind !=
size_type(-1))
7261 gmm::add(gmm::scaled(md.real_variable(vl[i], 0), half),
7262 gmm::scaled(md.real_variable(vl[i], 1), half),
7263 md.set_real_variable(vl[i], ind));
7264 md.set_default_iter_of_variable(vl[i], ind);
7269 for (
size_type i = 0; i < dl.size(); ++i) {
7270 bool is_uptodate = md.temporary_uptodate(dl[i], id_num, ind);
7271 if (!is_uptodate && ind !=
size_type(-1)) {
7272 gmm::add(gmm::scaled(md.real_variable(dl[i], 0), half),
7273 gmm::scaled(md.real_variable(dl[i], 1), half),
7274 md.set_real_variable(dl[i], ind));
7276 md.set_default_iter_of_variable(dl[i], ind);
7280 md.brick_call(ib, version, 0);
7281 if (pbr->is_linear()) {
7285 md.linear_brick_add_to_rhs(ib, 1, 1);
7288 md.reset_default_iter_of_variables(dl);
7289 if (!(pbr->is_linear()))
7290 md.reset_default_iter_of_variables(vl);
7293 virtual void asm_complex_tangent_terms
7294 (
const model &md,
size_type ib, model::complex_matlist &,
7295 std::vector<model::complex_veclist> &vectl,
7296 std::vector<model::complex_veclist> &vectl_sym,
7297 build_version version)
const {
7299 scalar_type half = scalar_type(1)/scalar_type(2);
7300 pbrick pbr = md.brick_pointer(ib);
7303 const model::varnamelist &vl = md.varnamelist_of_brick(ib);
7304 const model::varnamelist &dl = md.datanamelist_of_brick(ib);
7306 if (!(pbr->is_linear())) {
7307 for (
size_type i = 0; i < vl.size(); ++i) {
7308 bool is_uptodate = md.temporary_uptodate(vl[i], id_num, ind);
7309 if (!is_uptodate && ind !=
size_type(-1))
7310 gmm::add(gmm::scaled(md.complex_variable(vl[i], 0), half),
7311 gmm::scaled(md.complex_variable(vl[i], 1), half),
7312 md.set_complex_variable(vl[i], ind));
7313 md.set_default_iter_of_variable(vl[i], ind);
7318 for (
size_type i = 0; i < dl.size(); ++i) {
7319 bool is_uptodate = md.temporary_uptodate(dl[i], id_num, ind);
7320 if (!is_uptodate && ind !=
size_type(-1)) {
7321 gmm::add(gmm::scaled(md.complex_variable(dl[i], 0), half),
7322 gmm::scaled(md.complex_variable(dl[i], 1), half),
7323 md.set_complex_variable(dl[i], ind));
7325 md.set_default_iter_of_variable(dl[i], ind);
7329 md.brick_call(ib, version, 0);
7330 if (pbr->is_linear()) {
7334 md.linear_brick_add_to_rhs(ib, 1, 1);
7337 md.reset_default_iter_of_variables(dl);
7338 if (!(pbr->is_linear()))
7339 md.reset_default_iter_of_variables(vl);
7342 midpoint_dispatcher() : virtual_dispatcher(2)
7343 { id_num = act_counter(); }
7348 pdispatcher pdispatch = std::make_shared<midpoint_dispatcher>();
7349 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.