37 model::model(
bool comp_version) {
38 init(); complex_version = comp_version;
39 is_linear_ = is_symmetric_ = is_coercive_ =
true;
41 time_integration = 0; init_step =
false; time_step = scalar_type(1);
48 pstring s1 = std::make_shared<std::string>(
"Hess_u");
49 tree1.add_name(s1->c_str(), 6, 0, s1);
50 tree1.root->name =
"u";
51 tree1.root->op_type = GA_NAME;
52 tree1.root->node_type = GA_NODE_MACRO_PARAM;
54 tree1.root->nbc2 = ga_parse_prefix_operator(*s1);
55 tree1.root->nbc3 = ga_parse_prefix_test(*s1);
56 ga_macro gam1(
"Hess", tree1, 1);
57 macro_dict.add_macro(gam1);
60 pstring s2 = std::make_shared<std::string>(
"Div_u");
61 tree2.add_name(s2->c_str(), 5, 0, s2);
62 tree2.root->name =
"u";
63 tree2.root->op_type = GA_NAME;
64 tree2.root->node_type = GA_NODE_MACRO_PARAM;
66 tree2.root->nbc2 = ga_parse_prefix_operator(*s2);
67 tree2.root->nbc3 = ga_parse_prefix_test(*s2);
68 ga_macro gam2(
"Div", tree2, 1);
69 macro_dict.add_macro(gam2);
72 void model::var_description::set_size() {
74 v_num_var_iter.resize(n_iter);
75 v_num_iter.resize(n_iter);
76 size_type s = is_fem_dofs ? passociated_mf()->nb_dof()
77 : (imd ? imd->nb_filtered_index()
78 * imd->nb_tensor_elem()
83 complex_value[i].resize(s);
85 real_value[i].resize(s);
86 if (is_affine_dependent) {
88 affine_complex_value.resize(s);
90 affine_real_value.resize(s);
94 size_type model::var_description::add_temporary(gmm::uint64_type id_num) {
96 for (; nit < n_iter + n_temp_iter ; ++nit)
97 if (v_num_var_iter[nit] == id_num)
break;
98 if (nit >= n_iter + n_temp_iter) {
100 v_num_var_iter.resize(nit+1);
101 v_num_var_iter[nit] = id_num;
102 v_num_iter.resize(nit+1);
106 complex_value.resize(n_iter + n_temp_iter);
107 complex_value[nit].resize(s);
110 real_value.resize(n_iter + n_temp_iter);
111 real_value[nit].resize(s);
117 void model::var_description::clear_temporaries() {
121 complex_value.resize(n_iter);
123 real_value.resize(n_iter);
126 bool model::check_name_validity(
const std::string &name,
bool assert)
const {
128 if (variables.count(name) != 0) {
129 GMM_ASSERT1(!assert,
"Variable " << name <<
" already exists");
131 }
else if (variable_groups.count(name) != 0) {
133 name <<
" corresponds to an already existing group of "
138 name <<
" corresponds to an already existing macro");
140 }
else if (name.compare(
"X") == 0) {
141 GMM_ASSERT1(!assert,
"X is a reserved keyword of the generic "
142 "assembly language");
146 int ga_valid = ga_check_name_validity(name);
148 GMM_ASSERT1(!assert,
"Invalid variable name, corresponds to an "
149 "operator or function name of the generic assembly language");
151 }
else if (ga_valid == 2) {
152 GMM_ASSERT1(!assert,
"Invalid variable name having a reserved "
153 "prefix used by the generic assembly language");
155 }
else if (ga_valid == 3) {
156 std::string org_name = sup_previous_and_dot_to_varname(name);
157 if (org_name.size() < name.size() &&
158 variables.find(org_name) != variables.end()) {
160 "Dot and Previous are reserved prefix used for time "
161 "integration schemes");
166 bool valid = !name.empty() && isalpha(name[0]);
168 for (
size_type i = 1; i < name.size(); ++i)
169 if (!(isalnum(name[i]) || name[i] ==
'_')) valid =
false;
170 GMM_ASSERT1(!assert || valid,
171 "Illegal variable name : \"" << name <<
"\"");
176 std::string res_name = name;
177 bool valid = check_name_validity(res_name,
false);
178 for (
size_type i = 2; !valid && i < 50; ++i) {
180 m << name <<
'_' << i;
182 valid = check_name_validity(res_name,
false);
184 for (
size_type i = 2; !valid && i < 1000; ++i) {
186 m <<
"new_" << name <<
'_' << i;
188 valid = check_name_validity(res_name,
false);
190 GMM_ASSERT1(valid,
"Illegal variable name: " << name);
195 model::VAR_SET::const_iterator
196 model::find_variable(
const std::string &name)
const {
197 VAR_SET::const_iterator it = variables.find(name);
198 GMM_ASSERT1(it != variables.end(),
"Undefined variable " << name);
202 const model::var_description &
203 model::variable_description(
const std::string &name)
const {
204 return find_variable(name)->second;
207 std::string sup_previous_and_dot_to_varname(std::string v) {
208 if (!(v.compare(0, 8,
"Previous")) && (v[8] ==
'_' || v[9] ==
'_')) {
209 v = v.substr((v[8] ==
'_') ? 9 : 10);
211 if (!(v.compare(0, 3,
"Dot")) && (v[3] ==
'_' || v[4] ==
'_')) {
212 v = v.substr((v[3] ==
'_') ? 4 : 5);
219 return name.substr(0, PREFIX_OLD_LENGTH) ==
PREFIX_OLD;
223 return is_old(name) ? name.substr(PREFIX_OLD_LENGTH) : name;
227 if (
is_old(name))
return false;
228 VAR_SET::const_iterator it = find_variable(name);
229 if (!(it->second.is_variable))
return false;
230 if (it->second.is_affine_dependent)
231 it = variables.find(it->second.org_name);
232 return it->second.is_disabled;
236 if (
is_old(name))
return true;
237 VAR_SET::const_iterator it = find_variable(name);
238 if (it->second.is_affine_dependent)
239 it = variables.find(it->second.org_name);
240 return !(it->second.is_variable) || it->second.is_disabled;
244 return is_old(name) || !(variable_description(name).is_variable);
248 if (
is_old(name))
return false;
249 const auto &var_descr = variable_description(name);
250 return var_descr.is_internal && !var_descr.is_disabled;
253 bool model::is_affine_dependent_variable(
const std::string &name)
const {
254 return !(
is_old(name)) && variable_description(name).is_affine_dependent;
258 model::org_variable(
const std::string &name)
const {
259 GMM_ASSERT1(is_affine_dependent_variable(name),
260 "For affine dependent variables only");
261 return variable_description(name).org_name;
265 model::factor_of_variable(
const std::string &name)
const {
266 return variable_description(name).alpha;
269 void model::set_factor_of_variable(
const std::string &name, scalar_type a) {
270 VAR_SET::iterator it = variables.find(name);
271 GMM_ASSERT1(it != variables.end(),
"Undefined variable " << name);
272 if (it->second.alpha != a) {
273 it->second.alpha = a;
274 for (
auto &v_num : it->second.v_num_data) v_num = act_counter();
278 bool model::is_im_data(
const std::string &name)
const {
283 model::pim_data_of_variable(
const std::string &name)
const {
287 const gmm::uint64_type &
288 model::version_number_of_data_variable(
const std::string &name,
290 VAR_SET::const_iterator it = find_variable(name);
291 if (niter ==
size_type(-1)) niter = it->second.default_iter;
292 return it->second.v_num_data[niter];
297 if (act_size_to_be_done) actualize_sizes();
299 return gmm::vect_size(crhs);
300 else if (with_internal && gmm::vect_size(full_rrhs))
301 return gmm::vect_size(full_rrhs);
303 return gmm::vect_size(rrhs);
306 void model::resize_global_system()
const {
309 for (
auto &&v : variables)
310 if (v.second.is_variable) {
311 if (v.second.is_disabled)
312 v.second.I = gmm::sub_interval(0,0);
313 else if (!v.second.is_affine_dependent && !v.second.is_internal) {
314 v.second.I = gmm::sub_interval(full_size, v.second.size());
315 full_size += v.second.size();
320 for (
auto &&v : variables)
321 if (v.second.is_internal && !v.second.is_disabled) {
322 v.second.I = gmm::sub_interval(full_size, v.second.size());
323 full_size += v.second.size();
326 for (
auto &&v : variables)
327 if (v.second.is_affine_dependent) {
328 v.second.I = variables.find(v.second.org_name)->second.I;
332 if (complex_version) {
341 if (full_size > primary_size) {
343 gmm::resize(internal_rTM, full_size-primary_size, primary_size);
352 for (dal::bv_visitor ib(valid_bricks); !ib.finished(); ++ib)
353 for (
const term_description &term : bricks[ib].tlist)
354 if (term.is_global) {
355 bricks[ib].terms_to_be_computed =
true;
360 void model::actualize_sizes()
const {
362 bool actualized =
false;
363 getfem::local_guard lock = locks_.get_lock();
364 if (actualized)
return;
366 act_size_to_be_done =
false;
368 std::map<std::string, std::vector<std::string> > multipliers;
369 std::set<std::string> tobedone;
385 for (
auto &&v : variables) {
386 const std::string &vname = v.first;
387 var_description &vdescr = v.second;
388 if (vdescr.is_fem_dofs && !vdescr.is_affine_dependent) {
389 if ((vdescr.filter & VDESCRFILTER_CTERM)
390 || (vdescr.filter & VDESCRFILTER_INFSUP)) {
391 VAR_SET::iterator vfilt = variables.find(vdescr.filter_var);
392 GMM_ASSERT1(vfilt != variables.end(),
"The primal variable of the"
393 " multiplier does not exist : " << vdescr.filter_var);
394 GMM_ASSERT1(vfilt->second.is_fem_dofs,
"The primal variable of "
395 "the multiplier is not a fem variable");
396 multipliers[vdescr.filter_var].push_back(vname);
397 if (vdescr.v_num < vdescr.mf->version_number() ||
398 vdescr.v_num < vfilt->second.mf->version_number()) {
399 tobedone.insert(vdescr.filter_var);
402 switch (vdescr.filter) {
403 case VDESCRFILTER_NO:
404 if (vdescr.v_num < vdescr.mf->version_number()) {
406 vdescr.v_num = act_counter();
409 case VDESCRFILTER_REGION:
410 if (vdescr.v_num < vdescr.mf->version_number()) {
412 dor = vdescr.mf->dof_on_region(vdescr.filter_region);
413 vdescr.partial_mf->adapt(dor);
415 vdescr.v_num = act_counter();
423 && vdescr.v_num < vdescr.imd->version_number()) {
425 vdescr.v_num = act_counter();
429 for (
auto &&v : variables) {
430 var_description &vdescr = v.second;
431 if (vdescr.is_fem_dofs && !(vdescr.is_affine_dependent) &&
432 ((vdescr.filter & VDESCRFILTER_CTERM)
433 || (vdescr.filter & VDESCRFILTER_INFSUP))) {
434 if (tobedone.count(vdescr.filter_var)) {
439 dal::bit_vector alldof; alldof.add(0, vdescr.mf->nb_dof());
440 vdescr.partial_mf->adapt(alldof);
442 vdescr.v_num = act_counter();
447 resize_global_system();
449 for (
const std::string &vname : tobedone) {
456 const std::vector<std::string> &mults = multipliers[vname];
457 const var_description &vdescr = variable_description(vname);
459 gmm::col_matrix< gmm::rsvector<scalar_type> > MGLOB;
460 if (mults.size() > 1) {
462 for (
const std::string &mult : mults)
463 s += variable_description(mult).mf->nb_dof();
467 std::set<size_type> glob_columns;
468 for (
const std::string &multname : mults) {
469 var_description &multdescr = variables.find(multname)->second;
475 gmm::col_matrix< gmm::rsvector<scalar_type> >
476 MM(vdescr.associated_mf().nb_dof(), multdescr.mf->nb_dof());
477 bool termadded =
false;
479 if (multdescr.filter & VDESCRFILTER_CTERM) {
481 for (dal::bv_visitor ib(valid_bricks); !ib.finished(); ++ib) {
482 const brick_description &brick = bricks[ib];
484 bool cplx =
is_complex() && brick.pbr->is_complex();
486 if (brick.tlist.size() == 0) {
487 bool varc =
false, multc =
false;
488 for (
const std::string &var : brick.vlist) {
489 if (multname.compare(var) == 0) multc =
true;
490 if (vname.compare(var) == 0) varc =
true;
493 GMM_ASSERT1(!cplx,
"Sorry, not taken into account");
494 generic_expressions.clear();
495 brick.terms_to_be_computed =
true;
496 update_brick(ib, BUILD_MATRIX);
497 if (generic_expressions.size()) {
498 GMM_TRACE2(
"Generic assembly for actualize sizes");
501 accumulated_distro<decltype(rTM)> distro_rTM(rTM);
503 ga_workspace workspace(*
this);
504 for (
const auto &ge : generic_expressions)
505 workspace.add_expression(ge.expr, ge.mim, ge.region,
506 2, ge.secondary_domain);
507 workspace.set_assembled_matrix(distro_rTM);
508 workspace.assembly(2);
511 gmm::add(gmm::sub_matrix(rTM, vdescr.I, multdescr.I), MM);
512 gmm::add(gmm::transposed
513 (gmm::sub_matrix(rTM, multdescr.I, vdescr.I)), MM);
519 for (
size_type j = 0; j < brick.tlist.size(); ++j) {
520 const term_description &term = brick.tlist[j];
522 if (term.is_matrix_term) {
523 if (term.is_global) {
524 bool varc =
false, multc =
false;
525 for (
const std::string var : brick.vlist) {
526 if (multname.compare(var) == 0) multc =
true;
527 if (vname.compare(var) == 0) varc =
true;
530 GMM_ASSERT1(!cplx,
"Sorry, not taken into account");
531 generic_expressions.clear();
533 brick.terms_to_be_computed =
true;
534 update_brick(ib, BUILD_MATRIX);
537 gmm::add(gmm::sub_matrix(brick.rmatlist[j],
538 multdescr.I, vdescr.I),
540 gmm::add(gmm::transposed(gmm::sub_matrix
542 vdescr.I, multdescr.I)),
546 }
else if (multname.compare(term.var1) == 0 &&
547 vname.compare(term.var2) == 0) {
549 brick.terms_to_be_computed =
true;
550 update_brick(ib, BUILD_MATRIX);
554 gmm::add(gmm::transposed(gmm::real_part(brick.cmatlist[j])),
557 gmm::add(gmm::transposed(brick.rmatlist[j]), MM);
560 }
else if (multname.compare(term.var2) == 0 &&
561 vname.compare(term.var1) == 0) {
563 brick.terms_to_be_computed =
true;
564 update_brick(ib, BUILD_MATRIX);
568 gmm::add(gmm::real_part(brick.cmatlist[j]), MM);
570 gmm::add(brick.rmatlist[j], MM);
578 GMM_WARNING1(
"No term found to filter multiplier " << multname
579 <<
". Variable is cancelled");
580 }
else if (multdescr.filter & VDESCRFILTER_INFSUP) {
581 mesh_region rg(multdescr.filter_region);
582 multdescr.filter_mim->linked_mesh().intersect_with_mpi_region(rg);
584 *(multdescr.mf), rg);
587 MPI_SUM_SPARSE_MATRIX(MM);
592 std::set<size_type> columns;
594 if (columns.size() == 0)
595 GMM_WARNING1(
"Empty basis found for multiplier " << multname);
597 if (mults.size() > 1) {
598 gmm::copy(MM, gmm::sub_matrix
600 gmm::sub_interval(0, vdescr.associated_mf().nb_dof()),
601 gmm::sub_interval(s, multdescr.mf->nb_dof())));
603 glob_columns.insert(s + icol);
604 s += multdescr.mf->nb_dof();
606 dal::bit_vector kept;
609 if (multdescr.filter & VDESCRFILTER_REGION)
610 kept &= multdescr.mf->dof_on_region(multdescr.filter_region);
611 multdescr.partial_mf->adapt(kept);
612 multdescr.set_size();
613 multdescr.v_num = act_counter();
621 if (mults.size() > 1) {
629 for (
const std::string &multname : mults) {
630 var_description &multdescr = variables.find(multname)->second;
631 dal::bit_vector kept;
632 size_type nbdof = multdescr.mf->nb_dof();
633 for (
const size_type &icol : glob_columns)
634 if (icol >= s && icol < s + nbdof) kept.add(icol-s);
635 if (multdescr.filter & VDESCRFILTER_REGION)
636 kept &= multdescr.mf->dof_on_region(multdescr.filter_region);
637 multdescr.partial_mf->adapt(kept);
638 multdescr.set_size();
639 multdescr.v_num = act_counter();
640 s += multdescr.mf->nb_dof();
648 resize_global_system();
660 if (variables.size() == 0)
661 ost <<
"Model with no variable nor data" << endl;
663 ost <<
"List of model variables and data:" << endl;
664 for (
const auto &v : variables) {
665 const var_description &vdescr = v.second;
666 if (vdescr.is_variable) ost <<
"Variable ";
668 ost << std::setw(30) << std::left << v.first;
669 if (vdescr.n_iter == 1)
672 ost << std::setw(2) << std::right << vdescr.n_iter <<
" copies ";
673 if (vdescr.is_fem_dofs) ost <<
"fem dependant ";
674 else ost <<
"constant size ";
676 ost << std::setw(8) << std::right << si;
678 ost <<
" double" << ((si > 1) ?
"s." :
".");
679 if (vdescr.is_variable &&
682 if (vdescr.imd != 0) ost <<
"\t (is im_data)";
683 if (vdescr.is_affine_dependent) ost <<
"\t (is affine dependent)";
686 for (
const auto &vargroup : variable_groups) {
687 ost <<
"Variable group " << std::setw(30) << std::left
689 if (vargroup.second.size()) {
691 for (
const std::string vname : vargroup.second) {
692 if (!first) ost <<
",";
else first =
false;
696 }
else ost <<
" empty" << endl;
701 void model::listresiduals(std::ostream &ost)
const {
703 if (variables.size() == 0)
704 ost <<
"Model with no variable nor data" << endl;
707 for (
const auto &v : variables) {
708 if (v.second.is_variable) {
709 const model_real_plain_vector &rhs = v.second.is_internal
711 const gmm::sub_interval &II = interval_of_variable(v.first);
712 scalar_type res = gmm::vect_norm2(gmm::sub_vector(rhs, II));
713 if (!firstvar) cout <<
", ";
714 ost <<
"res_" << v.first <<
"= " << std::setw(11) << res;
724 bgeot::multi_index sizes(1);
730 const bgeot::multi_index &sizes,
732 check_name_validity(name);
733 variables.emplace(name, var_description(
true,
is_complex(), 0, 0, niter));
734 variables[name].qdims = sizes;
735 act_size_to_be_done =
true;
736 variables[name].set_size();
737 GMM_ASSERT1(variables[name].qdim(),
738 "Variables of null size are not allowed");
743 bgeot::multi_index sizes(1);
749 const bgeot::multi_index &sizes) {
750 GMM_ASSERT1(!(variables[name].is_fem_dofs),
751 "Cannot explicitly resize a fem variable or data");
752 GMM_ASSERT1(variables[name].imd == 0,
753 "Cannot explicitly resize an im data");
754 variables[name].qdims = sizes;
755 variables[name].set_size();
760 bgeot::multi_index sizes(1);
766 const bgeot::multi_index &sizes,
768 check_name_validity(name);
769 variables.emplace(name, var_description(
false,
is_complex(), 0, 0, niter));
770 variables[name].qdims = sizes;
771 GMM_ASSERT1(variables[name].qdim(),
"Data of null size are not allowed");
772 variables[name].set_size();
776 const base_matrix &M) {
779 GMM_ASSERT1(!(
is_complex()),
"Sorry, complex version to be done");
784 const base_complex_matrix &M) {
787 GMM_ASSERT1(!(
is_complex()),
"Sorry, complex version to be done");
792 const base_tensor &t) {
794 GMM_ASSERT1(!(
is_complex()),
"Sorry, complex version to be done");
799 const base_complex_tensor &t) {
801 GMM_ASSERT1(!(
is_complex()),
"Sorry, complex version to be done");
807 check_name_validity(name);
808 variables.emplace(name,
809 var_description(
true,
is_complex(), 0, &imd, niter));
810 variables[name].set_size();
812 act_size_to_be_done =
true;
818 variables[name].is_internal =
true;
823 check_name_validity(name);
824 variables.emplace(name,
825 var_description(
false,
is_complex(), 0, &imd, niter));
826 variables[name].set_size();
832 check_name_validity(name);
833 variables.emplace(name, var_description(
true,
is_complex(), &mf, 0, niter,
835 variables[name].set_size();
837 act_size_to_be_done =
true;
838 leading_dim = std::max(leading_dim, mf.
linked_mesh().dim());
844 check_name_validity(name);
845 variables.emplace(name, var_description(
true,
is_complex(), &mf, 0, niter,
846 VDESCRFILTER_REGION, region));
847 variables[name].set_size();
848 act_size_to_be_done =
true;
853 const std::string &org_name,
855 check_name_validity(name);
856 VAR_SET::const_iterator it = find_variable(org_name);
857 GMM_ASSERT1(it->second.is_variable && !(it->second.is_affine_dependent),
858 "The original variable should be a variable");
859 variables.emplace(name, variables[org_name]);
860 variables[name].is_affine_dependent =
true;
861 variables[name].org_name = org_name;
862 variables[name].alpha = alpha;
863 variables[name].set_size();
868 bgeot::multi_index sizes(1);
874 const bgeot::multi_index &sizes,
size_type niter) {
875 check_name_validity(name);
876 variables.emplace(name, var_description(
false,
is_complex(), &mf, 0, niter,
878 variables[name].qdims = sizes;
879 GMM_ASSERT1(variables[name].qdim(),
"Data of null size are not allowed");
880 variables[name].set_size();
885 const std::string &primal_name,
887 check_name_validity(name);
888 variables.emplace(name, var_description(
true,
is_complex(), &mf, 0, niter,
891 variables[name].set_size();
892 act_size_to_be_done =
true;
897 size_type region,
const std::string &primal_name,
899 check_name_validity(name);
900 variables.emplace(name, var_description(
true,
is_complex(), &mf, 0, niter,
901 VDESCRFILTER_REGION_CTERM, region,
903 variables[name].set_size();
904 act_size_to_be_done =
true;
909 const std::string &primal_name,
912 check_name_validity(name);
913 variables.emplace(name, var_description(
true,
is_complex(), &mf, 0, niter,
914 VDESCRFILTER_INFSUP, region,
916 variables[name].set_size();
917 act_size_to_be_done =
true;
927 VAR_SET::iterator it = variables.find(name);
928 GMM_ASSERT1(it != variables.end(),
"Undefined variable " << name);
929 it->second.is_disabled = !enabled;
930 for (
auto &&v : variables) {
931 if (((v.second.filter & VDESCRFILTER_INFSUP) ||
932 (v.second.filter & VDESCRFILTER_CTERM))
933 && name.compare(v.second.filter_var) == 0) {
934 v.second.is_disabled = !enabled;
936 if (v.second.is_variable && v.second.is_affine_dependent
937 && name.compare(v.second.org_name) == 0)
938 v.second.is_disabled = !enabled;
940 if (!act_size_to_be_done) resize_global_system();
948 check_name_validity(name.substr(0, name.find(
"(")));
949 macro_dict.add_macro(name, expr);
953 { macro_dict.del_macro(name); }
956 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
957 valid_bricks.del(ib);
958 active_bricks.del(ib);
960 for (
size_type i = 0; i < bricks[ib].mims.size(); ++i) {
961 const mesh_im *mim = bricks[ib].mims[i];
963 for (dal::bv_visitor ibb(valid_bricks); !ibb.finished(); ++ibb) {
964 for (
size_type j = 0; j < bricks[ibb].mims.size(); ++j)
965 if (bricks[ibb].mims[j] == mim) found =
true;
967 for (
const auto &v : variables) {
968 if (v.second.is_fem_dofs &&
969 (v.second.filter & VDESCRFILTER_INFSUP) &&
970 mim == v.second.filter_mim) found =
true;
972 if (!found) sup_dependency(*mim);
975 is_linear_ = is_symmetric_ = is_coercive_ =
true;
976 for (dal::bv_visitor ibb(valid_bricks); !ibb.finished(); ++ibb) {
977 is_linear_ = is_linear_ && bricks[ibb].pbr->is_linear();
978 is_symmetric_ = is_symmetric_ && bricks[ibb].pbr->is_symmetric();
979 is_coercive_ = is_coercive_ && bricks[ibb].pbr->is_coercive();
981 bricks[ib] = brick_description();
985 for (dal::bv_visitor ibb(valid_bricks); !ibb.finished(); ++ibb) {
986 for (
const auto &vname : bricks[ibb].vlist)
987 GMM_ASSERT1(varname.compare(vname),
988 "Cannot delete a variable which is still used by a brick");
989 for (
const auto &dname : bricks[ibb].dlist)
990 GMM_ASSERT1(varname.compare(dname),
991 "Cannot delete a data which is still used by a brick");
994 VAR_SET::const_iterator it = find_variable(varname);
996 if (it->second.is_fem_dofs) {
999 for(VAR_SET::iterator it2 = variables.begin();
1000 it2 != variables.end(); ++it2) {
1001 if (it != it2 && it2->second.is_fem_dofs && mf == it2->second.mf)
1004 if (!found) sup_dependency(*mf);
1006 if (it->second.filter & VDESCRFILTER_INFSUP) {
1007 const mesh_im *mim = it->second.filter_mim;
1009 for (dal::bv_visitor ibb(valid_bricks); !ibb.finished(); ++ibb) {
1010 for (
size_type j = 0; j < bricks[ibb].mims.size(); ++j)
1011 if (bricks[ibb].mims[j] == mim) found =
true;
1013 for (VAR_SET::iterator it2 = variables.begin();
1014 it2 != variables.end(); ++it2) {
1015 if (it != it2 && it2->second.is_fem_dofs &&
1016 (it2->second.filter & VDESCRFILTER_INFSUP) &&
1017 mim == it2->second.filter_mim) found =
true;
1019 if (!found) sup_dependency(*mim);
1023 if (it->second.imd != 0) sup_dependency(*(it->second.imd));
1025 variables.erase(varname);
1026 act_size_to_be_done =
true;
1030 const varnamelist &datanames,
1031 const termlist &terms,
1032 const mimlist &mims,
size_type region) {
1033 size_type ib = valid_bricks.first_false();
1035 for (
size_type i = 0; i < terms.size(); ++i)
1036 if (terms[i].is_global && terms[i].is_matrix_term && pbr->is_linear())
1037 GMM_ASSERT1(
false,
"Global linear matrix terms are not allowed");
1039 if (ib == bricks.size())
1040 bricks.push_back(brick_description(pbr, varnames, datanames, terms,
1043 bricks[ib] = brick_description(pbr, varnames, datanames, terms,
1045 active_bricks.add(ib);
1046 valid_bricks.add(ib);
1053 "Impossible to add a complex brick to a real model");
1055 bricks[ib].cmatlist.resize(terms.size());
1056 bricks[ib].cveclist[0].resize(terms.size());
1057 bricks[ib].cveclist_sym[0].resize(terms.size());
1059 bricks[ib].rmatlist.resize(terms.size());
1060 bricks[ib].rveclist[0].resize(terms.size());
1061 bricks[ib].rveclist_sym[0].resize(terms.size());
1063 is_linear_ = is_linear_ && pbr->is_linear();
1064 is_symmetric_ = is_symmetric_ && pbr->is_symmetric();
1065 is_coercive_ = is_coercive_ && pbr->is_coercive();
1067 for (
const auto &vname : varnames)
1068 GMM_ASSERT1(variables.count(vname),
1069 "Undefined model variable " << vname);
1071 for (
const auto &dname : datanames)
1072 GMM_ASSERT1(variables.count(dname),
1073 "Undefined model data or variable " << dname);
1079 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
1081 bricks[ib].mims.push_back(&mim);
1082 add_dependency(mim);
1086 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
1088 bricks[ib].tlist = terms;
1089 if (
is_complex() && bricks[ib].pbr->is_complex()) {
1090 bricks.back().cmatlist.resize(terms.size());
1091 bricks.back().cveclist[0].resize(terms.size());
1092 bricks.back().cveclist_sym[0].resize(terms.size());
1094 bricks.back().rmatlist.resize(terms.size());
1095 bricks.back().rveclist[0].resize(terms.size());
1096 bricks.back().rveclist_sym[0].resize(terms.size());
1101 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
1103 bricks[ib].vlist = vl;
1104 for (
const auto &v : vl)
1105 GMM_ASSERT1(variables.count(v),
"Undefined model variable " << v);
1109 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
1111 bricks[ib].dlist = dl;
1112 for (
const auto &v : dl)
1113 GMM_ASSERT1(variables.count(v),
"Undefined model variable " << v);
1117 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
1119 bricks[ib].mims = ml;
1120 for (
const auto &mim : ml) add_dependency(*mim);
1124 GMM_ASSERT1(valid_bricks[ib],
"Inexistent brick");
1126 bricks[ib].is_update_brick = flag;
1129 void model::set_time(scalar_type t,
bool to_init) {
1130 static const std::string varname(
"t");
1131 VAR_SET::iterator it = variables.find(varname);
1132 if (it == variables.end()) {
1135 GMM_ASSERT1(it->second.size() == 1,
"Time data should be of size 1");
1137 if (it == variables.end() || to_init) {
1145 scalar_type model::get_time() {
1146 static const std::string varname(
"t");
1147 set_time(scalar_type(0),
false);
1154 void model::call_init_affine_dependent_variables(
int version) {
1155 for (VAR_SET::iterator it = variables.begin();
1156 it != variables.end(); ++it) {
1157 var_description &vdescr = it->second;
1158 if (vdescr.is_variable && vdescr.ptsc) {
1160 vdescr.ptsc->init_affine_dependent_variables_precomputation(*
this);
1162 vdescr.ptsc->init_affine_dependent_variables(*
this);
1167 void model::shift_variables_for_time_integration() {
1168 for (VAR_SET::iterator it = variables.begin();
1169 it != variables.end(); ++it)
1170 if (it->second.is_variable && it->second.ptsc)
1171 it->second.ptsc->shift_variables(*
this);
1174 void model::add_time_integration_scheme(
const std::string &varname,
1175 ptime_scheme ptsc) {
1176 VAR_SET::iterator it = variables.find(varname);
1177 GMM_ASSERT1(it != variables.end(),
"Undefined variable " << varname);
1178 GMM_ASSERT1(it->second.is_variable && !(it->second.is_affine_dependent),
1179 "Cannot apply an integration scheme to " << varname);
1180 it->second.ptsc = ptsc;
1187 time_integration = 1;
1190 void model::copy_init_time_derivative() {
1192 for (VAR_SET::iterator it = variables.begin();
1193 it != variables.end(); ++it)
1194 if (it->second.is_variable && it->second.ptsc) {
1196 std::string name_v, name_previous_v;
1197 it->second.ptsc->time_derivative_to_be_initialized(name_v,
1200 if (name_v.size()) {
1218 class APIDECL first_order_theta_method_scheme
1219 :
public virtual_time_scheme {
1221 std::string U, U0, V, V0;
1226 virtual void init_affine_dependent_variables(model &md)
const {
1227 scalar_type dt = md.get_time_step();
1228 scalar_type a = scalar_type(1)/(theta*dt);
1229 scalar_type b = (scalar_type(1)-theta)/theta;
1230 md.set_factor_of_variable(V, a);
1231 if (md.is_complex()) {
1232 gmm::add(gmm::scaled(md.complex_variable(U0), -complex_type(a)),
1233 gmm::scaled(md.complex_variable(V0), -complex_type(b)),
1234 md.set_complex_constant_part(V));
1237 gmm::add(gmm::scaled(md.real_variable(U0), -a),
1238 gmm::scaled(md.real_variable(V0), -b),
1239 md.set_real_constant_part(V));
1244 virtual void init_affine_dependent_variables_precomputation(model &md)
1246 scalar_type dt = md.get_time_step();
1247 md.set_factor_of_variable(V, scalar_type(1)/dt);
1248 if (md.is_complex()) {
1249 gmm::copy(gmm::scaled(md.complex_variable(U0), -complex_type(1)/dt),
1250 md.set_complex_constant_part(V));
1253 gmm::copy(gmm::scaled(md.real_variable(U0), -scalar_type(1)/dt),
1254 md.set_real_constant_part(V));
1258 virtual void time_derivative_to_be_initialized
1259 (std::string &name_v, std::string &name_previous_v)
const
1260 {
if (theta != scalar_type(1)) { name_v = V; name_previous_v = V0; } }
1262 virtual void shift_variables(model &md)
const {
1263 if (md.is_complex()) {
1264 gmm::copy(md.complex_variable(U), md.set_complex_variable(U0));
1265 gmm::copy(md.complex_variable(V), md.set_complex_variable(V0));
1267 gmm::copy(md.real_variable(U), md.set_real_variable(U0));
1268 gmm::copy(md.real_variable(V), md.set_real_variable(V0));
1273 first_order_theta_method_scheme(model &md, std::string varname,
1276 U0 =
"Previous_" + U;
1278 V0 =
"Previous_Dot_" + U;
1280 GMM_ASSERT1(theta > scalar_type(0) && theta <= scalar_type(1),
1281 "Invalid value of theta parameter for the theta-method");
1283 if (!(md.variable_exists(V)))
1284 md.add_affine_dependent_variable(V, U);
1285 const mesh_fem *mf = md.pmesh_fem_of_variable(U);
1286 size_type s = md.is_complex() ? gmm::vect_size(md.complex_variable(U))
1287 : gmm::vect_size(md.real_variable(U));
1290 if (!(md.variable_exists(U0))) md.add_fem_data(U0, *mf);
1291 if (!(md.variable_exists(V0))) md.add_fem_data(V0, *mf);
1293 if (!(md.variable_exists(U0))) md.add_fixed_size_data(U0, s);
1294 if (!(md.variable_exists(V0))) md.add_fixed_size_data(V0, s);
1301 void add_theta_method_for_first_order(model &md,
const std::string &varname,
1302 scalar_type theta) {
1304 = std::make_shared<first_order_theta_method_scheme>(md, varname,theta);
1305 md.add_time_integration_scheme(varname, ptsc);
1314 class APIDECL second_order_theta_method_scheme
1315 :
public virtual_time_scheme {
1317 std::string U, U0, V, V0, A, A0;
1323 virtual void init_affine_dependent_variables(model &md)
const {
1324 scalar_type dt = md.get_time_step();
1325 md.set_factor_of_variable(V, scalar_type(1)/(theta*dt));
1326 md.set_factor_of_variable(A, scalar_type(1)/(theta*theta*dt*dt));
1327 if (md.is_complex()) {
1328 gmm::add(gmm::scaled(md.complex_variable(U0),
1329 -complex_type(1)/(theta*dt)),
1330 gmm::scaled(md.complex_variable(V0),
1331 -(complex_type(1)-complex_type(theta))/theta),
1332 md.set_complex_constant_part(V));
1333 gmm::add(gmm::scaled(md.complex_variable(U0),
1334 -complex_type(1)/(theta*theta*dt*dt)),
1335 gmm::scaled(md.complex_variable(A0),
1336 -(complex_type(1)-complex_type(theta))/theta),
1337 md.set_complex_constant_part(A));
1338 gmm::add(gmm::scaled(md.complex_variable(V0),
1339 -complex_type(1)/(theta*theta*dt)),
1340 md.set_complex_constant_part(A));
1344 gmm::add(gmm::scaled(md.real_variable(U0),
1345 -scalar_type(1)/(theta*dt)),
1346 gmm::scaled(md.real_variable(V0),
1347 -(scalar_type(1)-theta)/theta),
1348 md.set_real_constant_part(V));
1349 gmm::add(gmm::scaled(md.real_variable(U0),
1350 -scalar_type(1)/(theta*theta*dt*dt)),
1351 gmm::scaled(md.real_variable(A0),
1352 -(scalar_type(1)-theta)/theta),
1353 md.set_real_constant_part(A));
1354 gmm::add(gmm::scaled(md.real_variable(V0),
1355 -scalar_type(1)/(theta*theta*dt)),
1356 md.set_real_constant_part(A));
1363 virtual void init_affine_dependent_variables_precomputation(model &md)
1365 scalar_type dt = md.get_time_step();
1366 md.set_factor_of_variable(V, scalar_type(1)/dt);
1367 md.set_factor_of_variable(A, scalar_type(1)/(dt*dt));
1368 if (md.is_complex()) {
1369 gmm::copy(gmm::scaled(md.complex_variable(U0),
1370 -complex_type(1)/dt),
1371 md.set_complex_constant_part(V));
1372 gmm::add(gmm::scaled(md.complex_variable(U0),
1373 -complex_type(1)/(dt*dt)),
1374 gmm::scaled(md.complex_variable(V0),
1375 -complex_type(1)/dt),
1376 md.set_complex_constant_part(A));
1378 gmm::copy(gmm::scaled(md.real_variable(U0),
1379 -scalar_type(1)/dt),
1380 md.set_real_constant_part(V));
1381 gmm::add(gmm::scaled(md.real_variable(U0),
1382 -scalar_type(1)/(dt*dt)),
1383 gmm::scaled(md.real_variable(V0),
1384 -scalar_type(1)/dt),
1385 md.set_real_constant_part(A));
1389 virtual void time_derivative_to_be_initialized
1390 (std::string &name_v, std::string &name_previous_v)
const
1391 {
if (theta != scalar_type(1)) { name_v = A; name_previous_v = A0; } }
1393 virtual void shift_variables(model &md)
const {
1394 if (md.is_complex()) {
1395 gmm::copy(md.complex_variable(U), md.set_complex_variable(U0));
1396 gmm::copy(md.complex_variable(V), md.set_complex_variable(V0));
1397 gmm::copy(md.complex_variable(A), md.set_complex_variable(A0));
1399 gmm::copy(md.real_variable(U), md.set_real_variable(U0));
1400 gmm::copy(md.real_variable(V), md.set_real_variable(V0));
1401 gmm::copy(md.real_variable(A), md.set_real_variable(A0));
1406 second_order_theta_method_scheme(model &md, std::string varname,
1409 U0 =
"Previous_" + U;
1411 V0 =
"Previous_Dot_" + U;
1413 A0 =
"Previous_Dot2_" + U;
1415 GMM_ASSERT1(theta > scalar_type(0) && theta <= scalar_type(1),
1416 "Invalid value of theta parameter for the theta-method");
1418 if (!(md.variable_exists(V)))
1419 md.add_affine_dependent_variable(V, U);
1420 if (!(md.variable_exists(A)))
1421 md.add_affine_dependent_variable(A, U);
1423 const mesh_fem *mf = md.pmesh_fem_of_variable(U);
1424 size_type s = md.is_complex() ? gmm::vect_size(md.complex_variable(U))
1425 : gmm::vect_size(md.real_variable(U));
1428 if (!(md.variable_exists(U0))) md.add_fem_data(U0, *mf);
1429 if (!(md.variable_exists(V0))) md.add_fem_data(V0, *mf);
1430 if (!(md.variable_exists(A0))) md.add_fem_data(A0, *mf);
1432 if (!(md.variable_exists(U0))) md.add_fixed_size_data(U0, s);
1433 if (!(md.variable_exists(V0))) md.add_fixed_size_data(V0, s);
1434 if (!(md.variable_exists(A0))) md.add_fixed_size_data(A0, s);
1441 void add_theta_method_for_second_order(model &md,
const std::string &varname,
1442 scalar_type theta) {
1443 ptime_scheme ptsc = std::make_shared<second_order_theta_method_scheme>
1445 md.add_time_integration_scheme(varname, ptsc);
1455 class APIDECL Newmark_scheme
1456 :
public virtual_time_scheme {
1458 std::string U, U0, V, V0, A, A0;
1459 scalar_type beta, gamma;
1464 virtual void init_affine_dependent_variables(model &md)
const {
1465 scalar_type dt = md.get_time_step();
1466 scalar_type a0 = scalar_type(1)/(beta*dt*dt), a1 = dt*a0;
1467 scalar_type a2 = (scalar_type(1) - scalar_type(2)*beta)
1468 / (scalar_type(2)*beta);
1469 scalar_type b0 = gamma/(beta*dt), b1 = (beta-gamma)/beta;
1470 scalar_type b2 = dt*(scalar_type(1)-gamma/(scalar_type(2)*beta));
1472 md.set_factor_of_variable(V, b0);
1473 md.set_factor_of_variable(A, a0);
1474 if (md.is_complex()) {
1475 gmm::add(gmm::scaled(md.complex_variable(U0), -complex_type(b0)),
1476 gmm::scaled(md.complex_variable(V0), complex_type(b1)),
1477 md.set_complex_constant_part(V));
1478 gmm::add(gmm::scaled(md.complex_variable(A0), complex_type(b2)),
1479 md.set_complex_constant_part(V));
1480 gmm::add(gmm::scaled(md.complex_variable(U0), -complex_type(a0)),
1481 gmm::scaled(md.complex_variable(V0), -complex_type(a1)),
1482 md.set_complex_constant_part(A));
1483 gmm::add(gmm::scaled(md.complex_variable(A0), -complex_type(a2)),
1484 md.set_complex_constant_part(A));
1486 gmm::add(gmm::scaled(md.real_variable(U0), -b0),
1487 gmm::scaled(md.real_variable(V0), b1),
1488 md.set_real_constant_part(V));
1489 gmm::add(gmm::scaled(md.real_variable(A0), b2),
1490 md.set_real_constant_part(V));
1491 gmm::add(gmm::scaled(md.real_variable(U0), -a0),
1492 gmm::scaled(md.real_variable(V0), -a1),
1493 md.set_real_constant_part(A));
1494 gmm::add(gmm::scaled(md.real_variable(A0), -a2),
1495 md.set_real_constant_part(A));
1502 virtual void init_affine_dependent_variables_precomputation(model &md)
1504 scalar_type dt = md.get_time_step();
1505 md.set_factor_of_variable(V, scalar_type(1)/dt);
1506 md.set_factor_of_variable(A, scalar_type(1)/(dt*dt));
1507 if (md.is_complex()) {
1508 gmm::copy(gmm::scaled(md.complex_variable(U0),
1509 -complex_type(1)/dt),
1510 md.set_complex_constant_part(V));
1511 gmm::add(gmm::scaled(md.complex_variable(U0),
1512 -complex_type(1)/(dt*dt)),
1513 gmm::scaled(md.complex_variable(V0),
1514 -complex_type(1)/dt),
1515 md.set_complex_constant_part(A));
1517 gmm::copy(gmm::scaled(md.real_variable(U0),
1518 -scalar_type(1)/dt),
1519 md.set_real_constant_part(V));
1520 gmm::add(gmm::scaled(md.real_variable(U0),
1521 -scalar_type(1)/(dt*dt)),
1522 gmm::scaled(md.real_variable(V0),
1523 -scalar_type(1)/dt),
1524 md.set_real_constant_part(A));
1528 virtual void time_derivative_to_be_initialized
1529 (std::string &name_v, std::string &name_previous_v)
const {
1530 if (beta != scalar_type(0.5) || gamma != scalar_type(1))
1531 { name_v = A; name_previous_v = A0; }
1534 virtual void shift_variables(model &md)
const {
1535 if (md.is_complex()) {
1536 gmm::copy(md.complex_variable(U), md.set_complex_variable(U0));
1537 gmm::copy(md.complex_variable(V), md.set_complex_variable(V0));
1538 gmm::copy(md.complex_variable(A), md.set_complex_variable(A0));
1540 gmm::copy(md.real_variable(U), md.set_real_variable(U0));
1541 gmm::copy(md.real_variable(V), md.set_real_variable(V0));
1542 gmm::copy(md.real_variable(A), md.set_real_variable(A0));
1547 Newmark_scheme(model &md, std::string varname,
1548 scalar_type be, scalar_type ga) {
1550 U0 =
"Previous_" + U;
1552 V0 =
"Previous_Dot_" + U;
1554 A0 =
"Previous_Dot2_" + U;
1555 beta = be; gamma = ga;
1556 GMM_ASSERT1(beta > scalar_type(0) && beta <= scalar_type(1)
1557 && gamma >= scalar_type(0.5) && gamma <= scalar_type(1),
1558 "Invalid parameter values for the Newmark scheme");
1560 if (!(md.variable_exists(V)))
1561 md.add_affine_dependent_variable(V, U);
1562 if (!(md.variable_exists(A)))
1563 md.add_affine_dependent_variable(A, U);
1565 const mesh_fem *mf = md.pmesh_fem_of_variable(U);
1566 size_type s = md.is_complex() ? gmm::vect_size(md.complex_variable(U))
1567 : gmm::vect_size(md.real_variable(U));
1570 if (!(md.variable_exists(U0))) md.add_fem_data(U0, *mf);
1571 if (!(md.variable_exists(V0))) md.add_fem_data(V0, *mf);
1572 if (!(md.variable_exists(A0))) md.add_fem_data(A0, *mf);
1574 if (!(md.variable_exists(U0))) md.add_fixed_size_data(U0, s);
1575 if (!(md.variable_exists(V0))) md.add_fixed_size_data(V0, s);
1576 if (!(md.variable_exists(A0))) md.add_fixed_size_data(A0, s);
1583 void add_Newmark_scheme(model &md,
const std::string &varname,
1584 scalar_type beta, scalar_type gamma) {
1585 ptime_scheme ptsc = std::make_shared<Newmark_scheme>
1586 (md, varname, beta, gamma);
1587 md.add_time_integration_scheme(varname, ptsc);
1596 class APIDECL Houbolt_scheme
1597 :
public virtual_time_scheme {
1599 std::string U, U01, U02, U03, V, A;
1604 virtual void init_affine_dependent_variables(model &md)
const {
1605 scalar_type dt = md.get_time_step();
1606 scalar_type a0 = scalar_type(2)/(dt*dt);
1607 scalar_type a1 = scalar_type(5)/(dt*dt);
1608 scalar_type a2 = scalar_type(4)/(dt*dt);
1609 scalar_type a3 = scalar_type(1)/(dt*dt);
1610 scalar_type b0 = scalar_type(11)/(scalar_type(6)*dt);
1611 scalar_type b1 = scalar_type(18)/(scalar_type(6)*dt);
1612 scalar_type b2 = scalar_type(9)/(scalar_type(6)*dt);
1613 scalar_type b3 = scalar_type(2)/(scalar_type(6)*dt);
1615 md.set_factor_of_variable(V, b0);
1616 md.set_factor_of_variable(A, a0);
1617 if (md.is_complex()) {
1618 gmm::add(gmm::scaled(md.complex_variable(U01), -complex_type(b1)),
1619 gmm::scaled(md.complex_variable(U02), complex_type(b2)),
1620 md.set_complex_constant_part(V));
1621 gmm::add(gmm::scaled(md.complex_variable(U03), -complex_type(b3)),
1622 md.set_complex_constant_part(V));
1623 gmm::add(gmm::scaled(md.complex_variable(U01), -complex_type(a1)),
1624 gmm::scaled(md.complex_variable(U02), complex_type(a2)),
1625 md.set_complex_constant_part(A));
1626 gmm::add(gmm::scaled(md.complex_variable(U03), -complex_type(a3)),
1627 md.set_complex_constant_part(A));
1629 gmm::add(gmm::scaled(md.real_variable(U01), -b1),
1630 gmm::scaled(md.real_variable(U02), b2),
1631 md.set_real_constant_part(V));
1632 gmm::add(gmm::scaled(md.real_variable(U03), -b3),
1633 md.set_real_constant_part(V));
1634 gmm::add(gmm::scaled(md.real_variable(U01), -a1),
1635 gmm::scaled(md.real_variable(U02), a2),
1636 md.set_real_constant_part(A));
1637 gmm::add(gmm::scaled(md.real_variable(U03), -a3),
1638 md.set_real_constant_part(A));
1642 virtual void init_affine_dependent_variables_precomputation(model &md)
1647 virtual void time_derivative_to_be_initialized
1648 (std::string &name_v, std::string &name_previous_v)
const {
1650 (void) name_previous_v;
1653 virtual void shift_variables(model &md)
const {
1654 if (md.is_complex()) {
1655 gmm::copy(md.complex_variable(U02), md.set_complex_variable(U03));
1656 gmm::copy(md.complex_variable(U01), md.set_complex_variable(U02));
1657 gmm::copy(md.complex_variable(U), md.set_complex_variable(U01));
1659 gmm::copy(md.real_variable(U02), md.set_real_variable(U03));
1660 gmm::copy(md.real_variable(U01), md.set_real_variable(U02));
1661 gmm::copy(md.real_variable(U), md.set_real_variable(U01));
1666 Houbolt_scheme(model &md, std::string varname) {
1668 U01 =
"Previous_" + U;
1669 U02 =
"Previous2_" + U;
1670 U03 =
"Previous3_" + U;
1674 if (!(md.variable_exists(V)))
1675 md.add_affine_dependent_variable(V, U);
1676 if (!(md.variable_exists(A)))
1677 md.add_affine_dependent_variable(A, U);
1679 const mesh_fem *mf = md.pmesh_fem_of_variable(U);
1680 size_type s = md.is_complex() ? gmm::vect_size(md.complex_variable(U))
1681 : gmm::vect_size(md.real_variable(U));
1684 if (!(md.variable_exists(U01))) md.add_fem_data(U01, *mf);
1685 if (!(md.variable_exists(U02))) md.add_fem_data(U02, *mf);
1686 if (!(md.variable_exists(U03))) md.add_fem_data(U03, *mf);
1688 if (!(md.variable_exists(U01))) md.add_fixed_size_data(U01, s);
1689 if (!(md.variable_exists(U02))) md.add_fixed_size_data(U02, s);
1690 if (!(md.variable_exists(U03))) md.add_fixed_size_data(U03, s);
1697 void add_Houbolt_scheme(model &md,
const std::string &varname) {
1698 ptime_scheme ptsc = std::make_shared<Houbolt_scheme>
1700 md.add_time_integration_scheme(varname, ptsc);
1712 GMM_ASSERT1(valid_bricks[ibrick],
"Inexistent brick");
1713 pbrick pbr = bricks[ibrick].pbr;
1715 bricks[ibrick].pdispatch = pdispatch;
1718 = std::max(
size_type(1), pdispatch->nbrhs());
1720 gmm::resize(bricks[ibrick].coeffs, nbrhs);
1723 bricks[ibrick].cveclist.resize(nbrhs);
1724 bricks[ibrick].cveclist_sym.resize(nbrhs);
1726 bricks[ibrick].cveclist[k] = bricks[ibrick].cveclist[0];
1727 bricks[ibrick].cveclist_sym[k] = bricks[ibrick].cveclist_sym[0];
1730 bricks[ibrick].rveclist.resize(nbrhs);
1731 bricks[ibrick].rveclist_sym.resize(nbrhs);
1733 bricks[ibrick].rveclist[k] = bricks[ibrick].rveclist[0];
1734 bricks[ibrick].rveclist_sym[k] = bricks[ibrick].rveclist_sym[0];
1741 GMM_ASSERT1(valid_bricks[ind_brick],
"Inexistent brick");
1742 GMM_ASSERT1(ind_var < bricks[ind_brick].vlist.size(),
1743 "Inexistent brick variable");
1744 return bricks[ind_brick].vlist[ind_var];
1749 GMM_ASSERT1(valid_bricks[ind_brick],
"Inexistent brick");
1750 GMM_ASSERT1(ind_data < bricks[ind_brick].dlist.size(),
1751 "Inexistent brick data");
1752 return bricks[ind_brick].dlist[ind_data];
1756 if (valid_bricks.card() == 0)
1757 ost <<
"Model with no bricks" << endl;
1759 ost <<
"List of model bricks:" << endl;
1760 for (dal::bv_visitor i(valid_bricks); !i.finished(); ++i) {
1761 ost <<
"Brick " << std::setw(3) << std::right << i + base_id
1762 <<
" " << std::setw(20) << std::right
1763 << bricks[i].pbr->brick_name();
1764 if (!(active_bricks[i])) ost <<
" (deactivated)";
1765 if (bricks[i].pdispatch) ost <<
" (dispatched)";
1766 ost << endl <<
" concerned variables: " << bricks[i].vlist[0];
1767 for (
size_type j = 1; j < bricks[i].vlist.size(); ++j)
1768 ost <<
", " << bricks[i].vlist[j];
1770 ost <<
" brick with " << bricks[i].tlist.size() <<
" term";
1771 if (bricks[i].tlist.size() > 1) ost <<
"s";
1780 void model::brick_init(
size_type ib, build_version version,
1782 const brick_description &brick = bricks[ib];
1783 bool cplx =
is_complex() && brick.pbr->is_complex();
1786 for (
size_type j = 0; j < brick.tlist.size(); ++j) {
1787 const term_description &term = brick.tlist[j];
1788 bool isg = term.is_global;
1790 gmm::vect_size(crhs) : gmm::vect_size(rrhs);
1791 size_type nbd1 = isg ? nbgdof : variables[term.var1].size();
1792 size_type nbd2 = isg ? nbgdof : (term.is_matrix_term ?
1793 variables[term.var2].size() : 0);
1794 if (term.is_matrix_term &&
1795 (brick.pbr->is_linear() || (version | BUILD_MATRIX))) {
1796 if (version | BUILD_ON_DATA_CHANGE) {
1798 gmm::resize(brick.cmatlist[j], nbd1, nbd2);
1800 gmm::resize(brick.rmatlist[j], nbd1, nbd2);
1803 brick.cmatlist[j] = model_complex_sparse_matrix(nbd1, nbd2);
1805 brick.rmatlist[j] = model_real_sparse_matrix(nbd1, nbd2);
1808 if (brick.pbr->is_linear() || (version | BUILD_RHS)) {
1809 for (
size_type k = 0; k < brick.nbrhs; ++k) {
1811 if (k == rhs_ind)
gmm::clear(brick.cveclist[k][j]);
1813 if (term.is_symmetric && term.var1.compare(term.var2)) {
1814 if (k == rhs_ind)
gmm::clear(brick.cveclist_sym[k][j]);
1818 if (k == rhs_ind)
gmm::clear(brick.rveclist[k][j]);
1820 if (term.is_symmetric && term.var1.compare(term.var2)) {
1821 if (k == rhs_ind)
gmm::clear(brick.rveclist_sym[k][j]);
1830 void model::post_to_variables_step(){}
1832 void model::brick_call(
size_type ib, build_version version,
1835 const brick_description &brick = bricks[ib];
1836 bool cplx =
is_complex() && brick.pbr->is_complex();
1838 brick_init(ib, version, rhs_ind);
1842 brick.pbr->complex_pre_assembly_in_serial(*
this, ib, brick.vlist,
1843 brick.dlist, brick.mims,
1845 brick.cveclist[rhs_ind],
1846 brick.cveclist_sym[rhs_ind],
1847 brick.region, version);
1852 accumulated_distro<complex_matlist> cmatlist(brick.cmatlist);
1853 accumulated_distro<complex_veclist> cveclist(brick.cveclist[rhs_ind]);
1854 accumulated_distro<complex_veclist> cveclist_sym(brick.cveclist_sym[rhs_ind]);
1858 brick.pbr->asm_complex_tangent_terms(*
this, ib, brick.vlist,
1859 brick.dlist, brick.mims,
1863 brick.region, version);
1866 brick.pbr->complex_post_assembly_in_serial(*
this, ib, brick.vlist,
1867 brick.dlist, brick.mims,
1869 brick.cveclist[rhs_ind],
1870 brick.cveclist_sym[rhs_ind],
1871 brick.region, version);
1873 if (brick.is_update_brick)
1875 for (
auto &&mat : brick.cmatlist)
1878 for (
auto &&vecs : brick.cveclist)
1879 for (
auto &&vec : vecs)
1882 for (
auto &&vecs : brick.cveclist_sym)
1883 for (
auto &&vec : vecs)
1889 brick.pbr->real_pre_assembly_in_serial(*
this, ib, brick.vlist,
1890 brick.dlist, brick.mims,
1892 brick.rveclist[rhs_ind],
1893 brick.rveclist_sym[rhs_ind],
1894 brick.region, version);
1897 accumulated_distro<real_matlist> rmatlist(brick.rmatlist);
1898 accumulated_distro<real_veclist> rveclist(brick.rveclist[rhs_ind]);
1899 accumulated_distro<real_veclist> rveclist_sym(brick.rveclist_sym[rhs_ind]);
1903 brick.pbr->asm_real_tangent_terms(*
this, ib, brick.vlist,
1904 brick.dlist, brick.mims,
1912 brick.pbr->real_post_assembly_in_serial(*
this, ib, brick.vlist,
1913 brick.dlist, brick.mims,
1915 brick.rveclist[rhs_ind],
1916 brick.rveclist_sym[rhs_ind],
1917 brick.region, version);
1919 if (brick.is_update_brick)
1921 for (
auto &&mat : brick.rmatlist)
1924 for (
auto &&vecs : brick.rveclist)
1925 for (
auto &&vec : vecs)
1928 for (
auto &&vecs : brick.rveclist_sym)
1929 for (
auto &&vec : vecs)
1936 void model::set_dispatch_coeff() {
1937 for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) {
1938 brick_description &brick = bricks[ib];
1939 if (brick.pdispatch)
1940 brick.pdispatch->set_dispatch_coeff(*
this, ib);
1946 context_check();
if (act_size_to_be_done) actualize_sizes();
1947 for (
auto && v : variables) v.second.clear_temporaries();
1949 set_dispatch_coeff();
1951 for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) {
1952 brick_description &brick = bricks[ib];
1953 if (brick.pdispatch) {
1955 brick.pdispatch->next_complex_iter(*
this, ib, brick.vlist,
1957 brick.cmatlist, brick.cveclist,
1958 brick.cveclist_sym,
true);
1960 brick.pdispatch->next_real_iter(*
this, ib, brick.vlist, brick.dlist,
1961 brick.rmatlist, brick.rveclist,
1962 brick.rveclist_sym,
true);
1968 context_check();
if (act_size_to_be_done) actualize_sizes();
1969 set_dispatch_coeff();
1971 for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) {
1972 brick_description &brick = bricks[ib];
1973 if (brick.pdispatch) {
1975 brick.pdispatch->next_complex_iter(*
this, ib, brick.vlist,
1977 brick.cmatlist, brick.cveclist,
1978 brick.cveclist_sym,
false);
1980 brick.pdispatch->next_real_iter(*
this, ib, brick.vlist, brick.dlist,
1981 brick.rmatlist, brick.rveclist,
1982 brick.rveclist_sym,
false);
1986 for (
auto &&v : variables)
1987 for (
size_type i = 1; i < v.second.n_iter; ++i) {
1989 gmm::copy(v.second.complex_value[i-1], v.second.complex_value[i]);
1991 gmm::copy(v.second.real_value[i-1], v.second.real_value[i]);
1992 v.second.v_num_data[i] = act_counter();
1996 bool model::is_var_newer_than_brick(
const std::string &varname,
1998 const brick_description &brick = bricks[ib];
1999 var_description &vd = variables[varname];
2000 if (niter ==
size_type(-1)) niter = vd.default_iter;
2001 return (vd.v_num > brick.v_num || vd.v_num_data[niter] > brick.v_num);
2004 bool model::is_var_mf_newer_than_brick(
const std::string &varname,
2006 const brick_description &brick = bricks[ib];
2007 var_description &vd = variables[varname];
2008 return (vd.v_num > brick.v_num);
2011 bool model::is_mim_newer_than_brick(
const mesh_im &im,
2013 const brick_description &brick = bricks[ib];
2014 return (im.version_number() > brick.v_num);
2017 void model::define_variable_group(
const std::string &group_name,
2018 const std::vector<std::string> &nl) {
2020 "variables cannot be the same as a variable name");
2022 std::set<const mesh *> ms;
2023 bool is_data_ =
false;
2024 for (
size_type i = 0; i < nl.size(); ++i) {
2029 "It is not possible to mix variables and data in a group");
2032 "All variables in a group have to exist in the model");
2034 GMM_ASSERT1(mf,
"Variables in a group should be fem variables");
2035 GMM_ASSERT1(ms.find(&(mf->linked_mesh())) == ms.end(),
2036 "Two variables in a group cannot share the same mesh");
2037 ms.insert(&(mf->linked_mesh()));
2039 variable_groups[group_name] = nl;
2042 void model::add_assembly_assignments(
const std::string &varname,
2045 GMM_ASSERT1(order < 3 || order ==
size_type(-1),
"Bad order value");
2046 const im_data *imd = pim_data_of_variable(varname);
2047 GMM_ASSERT1(imd != 0,
"Only applicable to im_data");
2048 assignement_desc as;
2049 as.varname = varname; as.expr = expr; as.region = rg; as.order = order;
2051 assignments.push_back(as);
2054 void model::add_temporaries(
const varnamelist &vl,
2055 gmm::uint64_type id_num)
const {
2056 for (
size_type i = 0; i < vl.size(); ++i) {
2057 var_description &vd = variables[vl[i]];
2058 if (vd.n_iter > 1) {
2059 vd.add_temporary(id_num);
2064 bool model::temporary_uptodate(
const std::string &varname,
2065 gmm::uint64_type id_num,
2067 var_description &vd = variables[varname];
2069 for (; ind < vd.n_iter + vd.n_temp_iter ; ++ind) {
2070 if (vd.v_num_var_iter[ind] == id_num)
break;
2072 if (ind < vd.n_iter + vd.n_temp_iter) {
2073 if (vd.v_num_iter[ind] <= vd.v_num_data[vd.default_iter]) {
2074 vd.v_num_iter[ind] = act_counter();
2083 void model::set_default_iter_of_variable(
const std::string &varname,
2086 var_description &vd = variables[varname];
2087 GMM_ASSERT1(ind < vd.n_iter + vd.n_temp_iter,
2088 "Inexistent iteration " << ind);
2089 vd.default_iter = ind;
2093 void model::reset_default_iter_of_variables(
const varnamelist &vl)
const {
2094 for (
size_type i = 0; i < vl.size(); ++i)
2095 variables[vl[i]].default_iter = 0;
2098 const model_real_sparse_matrix &
2100 GMM_ASSERT1(bricks[ib].tlist[iterm].is_matrix_term,
2101 "Not a matrix term !");
2102 GMM_ASSERT1(bricks[ib].pbr->is_linear(),
"Nonlinear term !");
2103 return bricks[ib].rmatlist[iterm];
2106 const model_complex_sparse_matrix &
2108 GMM_ASSERT1(bricks[ib].tlist[iterm].is_matrix_term,
2109 "Not a matrix term !");
2110 GMM_ASSERT1(bricks[ib].pbr->is_linear(),
"Nonlinear term !");
2111 return bricks[ib].cmatlist[iterm];
2115 void model::update_brick(
size_type ib, build_version version)
const {
2116 const brick_description &brick = bricks[ib];
2117 bool cplx =
is_complex() && brick.pbr->is_complex();
2118 bool tobecomputed = brick.terms_to_be_computed
2119 || brick.pbr->is_to_be_computed_each_time()
2120 || !(brick.pbr->is_linear());
2123 if (!tobecomputed ) {
2124 for (
size_type i = 0; i < brick.vlist.size() && !tobecomputed; ++i) {
2125 var_description &vd = variables[brick.vlist[i]];
2126 if (vd.v_num > brick.v_num) tobecomputed =
true;
2131 for (
size_type i = 0; i < brick.dlist.size() && !tobecomputed; ++i) {
2132 var_description &vd = variables[brick.dlist[i]];
2133 if (vd.v_num > brick.v_num || vd.v_num_data[vd.default_iter] > brick.v_num) {
2134 tobecomputed =
true;
2135 version = build_version(version | BUILD_ON_DATA_CHANGE);
2140 if (!tobecomputed ) {
2141 for (
size_type i = 0; i < brick.mims.size() && !tobecomputed; ++i) {
2142 if (brick.mims[i]->version_number() > brick.v_num) tobecomputed =
true;
2147 brick.external_load = scalar_type(0);
2149 if (!(brick.pdispatch))
2150 { brick_call(ib, version, 0); }
2153 brick.pdispatch->asm_complex_tangent_terms
2154 (*
this, ib, brick.cmatlist, brick.cveclist, brick.cveclist_sym,
2157 brick.pdispatch->asm_real_tangent_terms
2158 (*
this, ib, brick.rmatlist, brick.rveclist, brick.rveclist_sym,
2161 brick.v_num = act_counter();
2164 if (brick.pbr->is_linear()) brick.terms_to_be_computed =
false;
2171 const brick_description &brick = bricks[ib];
2172 if (brick.pbr->is_linear()) {
2174 bool cplx =
is_complex() && brick.pbr->is_complex();
2176 for (
size_type j = 0; j < brick.tlist.size(); ++j) {
2177 const term_description &term = brick.tlist[j];
2178 bool isg = term.is_global;
2181 size_type n_iter_1 = n_iter, n_iter_2 = n_iter;
2183 n_iter_1 = variables[term.var1].default_iter;
2184 if (term.is_matrix_term)
2185 n_iter_2 = variables[term.var2].default_iter;
2190 if (term.is_matrix_term) {
2193 model_complex_plain_vector V(nbgdof);
2194 for (VAR_SET::iterator it = variables.begin();
2195 it != variables.end(); ++it)
2196 if (it->second.is_variable) {
2198 ? it->second.default_iter : n_iter;
2199 gmm::copy(it->second.complex_value[n_iter_i],
2200 gmm::sub_vector(V, it->second.I));
2204 gmm::scaled(V, complex_type(-1)),
2205 brick.cveclist[ind_data][j]);
2209 gmm::scaled(variables[term.var2].complex_value[n_iter_2],
2211 brick.cveclist[ind_data][j]);
2215 model_real_plain_vector V(nbgdof);
2216 for (VAR_SET::iterator it = variables.begin();
2217 it != variables.end(); ++it)
2218 if (it->second.is_variable) {
2220 ? it->second.default_iter : n_iter;
2221 gmm::copy(it->second.real_value[n_iter_i],
2222 gmm::sub_vector(V, it->second.I));
2225 (brick.rmatlist[j], gmm::scaled(V, scalar_type(-1)),
2226 brick.rveclist[ind_data][j]);
2230 gmm::scaled(variables[term.var2].real_value[n_iter_2],
2231 scalar_type(-1)), brick.rveclist[ind_data][j]);
2234 if (term.is_symmetric && term.var1.compare(term.var2)) {
2237 (gmm::conjugated(brick.cmatlist[j]),
2238 gmm::scaled(variables[term.var1].complex_value[n_iter_1],
2240 brick.cveclist_sym[ind_data][j]);
2243 (gmm::transposed(brick.rmatlist[j]),
2244 gmm::scaled(variables[term.var1].real_value[n_iter_1],
2246 brick.rveclist_sym[ind_data][j]);
2253 void model::update_affine_dependent_variables() {
2254 for (VAR_SET::iterator it = variables.begin(); it != variables.end(); ++it)
2255 if (it->second.is_affine_dependent) {
2256 VAR_SET::iterator it2 = variables.find(it->second.org_name);
2257 if (it->second.size() != it2->second.size())
2258 it->second.set_size();
2259 if (it->second.is_complex) {
2260 gmm::add(gmm::scaled(it2->second.complex_value[0],
2261 complex_type(it->second.alpha)),
2262 it->second.affine_complex_value,
2263 it->second.complex_value[0]);
2265 gmm::add(gmm::scaled(it2->second.real_value[0], it->second.alpha),
2266 it->second.affine_real_value, it->second.real_value[0]);
2268 it->second.v_num = std::max(it->second.v_num, it2->second.v_num);
2269 for (
size_type i = 0; i < it->second.n_iter; ++i)
2271 it->second.v_num_data[i] = std::max(it->second.v_num_data[i],
2272 it2->second.v_num_data[i]);
2283 GMM_ASSERT1(mf,
"Works only with fem variables.");
2287 for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) {
2288 brick_description &brick = bricks[ib];
2290 bool detected =
false;
2291 for (
size_type i = 0; i < brick.vlist.size(); ++i)
2292 if (brick.vlist[i].compare(varname) == 0)
2293 { detected =
true;
break; }
2295 if (detected && brick.mims.size()) {
2297 for (
auto &pmim : brick.mims)
2300 pmim->linked_mesh()));
2301 GMM_ASSERT1(ifo >= 0,
2302 "The given region is only partially covered by "
2303 "region of brick \"" << brick.pbr->brick_name()
2304 <<
"\". Please subdivise the region");
2306 std::string expr = brick.pbr->declare_volume_assembly_string
2307 (*
this, ib, brick.vlist, brick.dlist);
2309 ga_workspace workspace(*
this);
2310 size_type order = workspace.add_expression
2311 (expr, dummy_mim, region);
2312 GMM_ASSERT1(order <= 1,
"Wrong order for a Neumann term");
2313 expr = workspace.extract_Neumann_term(varname);
2316 result +=
" + " + workspace.extract_Neumann_term(varname);
2318 result = workspace.extract_Neumann_term(varname);
2330 GMM_ASSERT1(version != BUILD_ON_DATA_CHANGE,
2331 "Invalid assembly version BUILD_ON_DATA_CHANGE");
2332 GMM_ASSERT1(version != BUILD_WITH_LIN,
2333 "Invalid assembly version BUILD_WITH_LIN");
2334 GMM_ASSERT1(version != BUILD_WITH_INTERNAL,
2335 "Invalid assembly version BUILD_WITH_INTERNAL");
2336 #if GETFEM_PARA_LEVEL > 1
2337 double t_ref = MPI_Wtime();
2340 context_check();
if (act_size_to_be_done) actualize_sizes();
2342 if (version & BUILD_MATRIX) gmm::clear(cTM);
2343 if (version & BUILD_RHS) gmm::clear(crhs);
2346 if (version & BUILD_MATRIX) gmm::clear(rTM);
2347 if (version & BUILD_RHS) gmm::clear(rrhs);
2349 clear_dof_constraints();
2350 generic_expressions.clear();
2351 update_affine_dependent_variables();
2353 if (version & BUILD_RHS) approx_external_load_ = scalar_type(0);
2355 for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) {
2357 brick_description &brick = bricks[ib];
2360 bool auto_disabled_brick =
true;
2361 for (
size_type j = 0; j < brick.vlist.size(); ++j) {
2363 auto_disabled_brick =
false;
2365 if (auto_disabled_brick)
continue;
2367 update_brick(ib, version);
2369 bool cplx =
is_complex() && brick.pbr->is_complex();
2371 scalar_type coeff0 = scalar_type(1);
2372 if (brick.pdispatch) coeff0 = brick.matrix_coeff;
2376 for (
size_type j = 0; j < brick.tlist.size(); ++j) {
2377 term_description &term = brick.tlist[j];
2378 bool isg = term.is_global, isprevious =
false;
2380 scalar_type alpha = coeff0, alpha1 = coeff0, alpha2 = coeff0;
2381 gmm::sub_interval I1(0,nbgdof), I2(0,nbgdof);
2382 var_description *var1=
nullptr, *var2=
nullptr;
2384 VAR_SET::iterator it1 = variables.find(term.var1);
2385 GMM_ASSERT1(it1 != variables.end(),
"Internal error");
2386 var1 = &(it1->second);
2387 GMM_ASSERT1(var1->is_variable,
"Assembly of data not allowed");
2389 if (term.is_matrix_term) {
2390 VAR_SET::iterator it2 = variables.find(term.var2);
2391 GMM_ASSERT1(it2 != variables.end(),
"Internal error");
2392 var2 = &(it2->second);
2394 if (!(var2->is_variable)) {
2395 std::string vorgname = sup_previous_and_dot_to_varname(term.var2);
2396 VAR_SET::iterator it3 = variables.find(vorgname);
2397 GMM_ASSERT1(it3->second.is_variable,
2398 "Assembly of data not allowed");
2402 alpha *= var1->alpha * var2->alpha;
2403 alpha1 *= var1->alpha;
2404 alpha2 *= var2->alpha;
2409 if (term.is_matrix_term && (version & BUILD_MATRIX) && !isprevious
2410 && (isg || (!(var1->is_disabled) && !(var2->is_disabled)))) {
2411 gmm::add(gmm::scaled(brick.cmatlist[j], alpha),
2412 gmm::sub_matrix(cTM, I1, I2));
2413 if (term.is_symmetric && I1.first() != I2.first()) {
2414 gmm::add(gmm::scaled(gmm::transposed(brick.cmatlist[j]), alpha),
2415 gmm::sub_matrix(cTM, I2, I1));
2418 if (version & BUILD_RHS) {
2419 if (isg || !(var1->is_disabled)) {
2420 if (brick.pdispatch) {
2421 for (
size_type k = 0; k < brick.nbrhs; ++k)
2422 gmm::add(gmm::scaled(brick.cveclist[k][j],
2424 gmm::sub_vector(crhs, I1));
2427 gmm::add(gmm::scaled(brick.cveclist[0][j],
2428 complex_type(alpha1)),
2429 gmm::sub_vector(crhs, I1));
2432 if (term.is_matrix_term && brick.pbr->is_linear() &&
is_linear()) {
2433 if (var2->is_affine_dependent
2434 && !(var1->is_disabled))
2435 gmm::mult_add(brick.cmatlist[j],
2436 gmm::scaled(var2->affine_complex_value,
2437 complex_type(-alpha1)),
2438 gmm::sub_vector(crhs, I1));
2439 if (term.is_symmetric && I1.first() != I2.first()
2440 && var1->is_affine_dependent
2441 && !(var2->is_disabled)) {
2442 gmm::mult_add(gmm::conjugated(brick.cmatlist[j]),
2443 gmm::scaled(var1->affine_complex_value,
2444 complex_type(-alpha2)),
2445 gmm::sub_vector(crhs, I2));
2448 if (term.is_matrix_term && brick.pbr->is_linear()
2449 && (!
is_linear() || (version & BUILD_WITH_LIN))) {
2450 if (!(var1->is_disabled))
2451 gmm::mult_add(brick.cmatlist[j],
2452 gmm::scaled(var2->complex_value[0],
2453 complex_type(-alpha1)),
2454 gmm::sub_vector(crhs, I1));
2456 if (term.is_symmetric && I1.first() != I2.first() &&
2457 !(var2->is_disabled)) {
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));
2468 if (brick.pbr->is_linear()
2469 && (!
is_linear() || (version & BUILD_WITH_LIN))) {
2470 gmm::mult_add(gmm::conjugated(brick.cmatlist[j]),
2471 gmm::scaled(var1->complex_value[0],
2472 complex_type(-alpha2)),
2473 gmm::sub_vector(crhs, I2));
2478 if (term.is_matrix_term && (version & BUILD_MATRIX) && !isprevious
2479 && (isg || (!(var1->is_disabled) && !(var2->is_disabled)))) {
2480 gmm::add(gmm::scaled(brick.rmatlist[j], alpha),
2481 gmm::sub_matrix(cTM, I1, I2));
2482 if (term.is_symmetric && I1.first() != I2.first()) {
2483 gmm::add(gmm::scaled(gmm::transposed(brick.rmatlist[j]), alpha),
2484 gmm::sub_matrix(cTM, I2, I1));
2487 if (version & BUILD_RHS) {
2488 if (isg || !(var1->is_disabled)) {
2489 if (brick.pdispatch) {
2490 for (
size_type k = 0; k < brick.nbrhs; ++k)
2491 gmm::add(gmm::scaled(brick.rveclist[k][j],
2493 gmm::sub_vector(crhs, I1));
2496 gmm::add(gmm::scaled(brick.rveclist[0][j], alpha1),
2497 gmm::sub_vector(crhs, I1));
2500 if (term.is_matrix_term && brick.pbr->is_linear() &&
is_linear()) {
2501 if (var2->is_affine_dependent
2502 && !(var1->is_disabled))
2503 gmm::mult_add(brick.rmatlist[j],
2504 gmm::scaled(var2->affine_complex_value,
2505 complex_type(-alpha1)),
2506 gmm::sub_vector(crhs, I1));
2507 if (term.is_symmetric && I1.first() != I2.first()
2508 && var1->is_affine_dependent
2509 && !(var2->is_disabled)) {
2510 gmm::mult_add(gmm::transposed(brick.rmatlist[j]),
2511 gmm::scaled(var1->affine_complex_value,
2512 complex_type(-alpha2)),
2513 gmm::sub_vector(crhs, I2));
2516 if (term.is_matrix_term && brick.pbr->is_linear()
2517 && (!
is_linear() || (version & BUILD_WITH_LIN))) {
2518 if (!(var1->is_disabled))
2519 gmm::mult_add(brick.rmatlist[j],
2520 gmm::scaled(var2->complex_value[0],
2521 complex_type(-alpha1)),
2522 gmm::sub_vector(crhs, I1));
2524 if (term.is_symmetric && I1.first() != I2.first() &&
2525 !(var2->is_disabled)) {
2526 if (brick.pdispatch) {
2527 for (
size_type k = 0; k < brick.nbrhs; ++k)
2528 gmm::add(gmm::scaled(brick.rveclist_sym[k][j],
2530 gmm::sub_vector(crhs, I2));
2533 gmm::add(gmm::scaled(brick.rveclist_sym[0][j], alpha2),
2534 gmm::sub_vector(crhs, I2));
2536 if (brick.pbr->is_linear()
2537 && (!
is_linear() || (version & BUILD_WITH_LIN))) {
2538 gmm::mult_add(gmm::transposed(brick.rmatlist[j]),
2539 gmm::scaled(var1->complex_value[0],
2540 complex_type(-alpha2)),
2541 gmm::sub_vector(crhs, I2));
2546 if (term.is_matrix_term && (version & BUILD_MATRIX) && !isprevious
2547 && (isg || (!(var1->is_disabled)
2548 && !(var2->is_disabled)))) {
2549 gmm::add(gmm::scaled(brick.rmatlist[j], alpha),
2550 gmm::sub_matrix(rTM, I1, I2));
2551 if (term.is_symmetric && I1.first() != I2.first()) {
2552 gmm::add(gmm::scaled(gmm::transposed(brick.rmatlist[j]), alpha),
2553 gmm::sub_matrix(rTM, I2, I1));
2556 if (version & BUILD_RHS) {
2558 if (isg || !(var1->is_disabled)) {
2559 if (brick.pdispatch) {
2560 for (
size_type k = 0; k < brick.nbrhs; ++k)
2561 gmm::add(gmm::scaled(brick.rveclist[k][j],
2563 gmm::sub_vector(rrhs, I1));
2566 gmm::add(gmm::scaled(brick.rveclist[0][j], alpha1),
2567 gmm::sub_vector(rrhs, I1));
2569 if (!(var1->is_disabled)) {
2571 if (term.is_matrix_term && brick.pbr->is_linear() &&
is_linear()
2572 && var2->is_affine_dependent)
2573 gmm::mult_add(brick.rmatlist[j],
2574 gmm::scaled(var2->affine_real_value, -alpha1),
2575 gmm::sub_vector(rrhs, I1));
2577 if (term.is_matrix_term && brick.pbr->is_linear()
2578 && (!
is_linear() || (version & BUILD_WITH_LIN)))
2579 gmm::mult_add(brick.rmatlist[j],
2580 gmm::scaled(var2->real_value[0], -alpha1),
2581 gmm::sub_vector(rrhs, I1));
2584 if (term.is_symmetric && I1.first() != I2.first() &&
2585 !(var2->is_disabled)) {
2586 if (brick.pdispatch) {
2587 for (
size_type k = 0; k < brick.nbrhs; ++k)
2588 gmm::add(gmm::scaled(brick.rveclist_sym[k][j],
2590 gmm::sub_vector(rrhs, I2));
2593 gmm::add(gmm::scaled(brick.rveclist_sym[0][j], alpha2),
2594 gmm::sub_vector(rrhs, I2));
2596 if (term.is_matrix_term && brick.pbr->is_linear() &&
is_linear()
2597 && var1->is_affine_dependent)
2598 gmm::mult_add(gmm::transposed(brick.rmatlist[j]),
2599 gmm::scaled(var1->affine_real_value, -alpha2),
2600 gmm::sub_vector(rrhs, I2));
2602 if (brick.pbr->is_linear()
2603 && (!
is_linear() || (version & BUILD_WITH_LIN)))
2604 gmm::mult_add(gmm::transposed(brick.rmatlist[j]),
2605 gmm::scaled(var1->real_value[0], -alpha2),
2606 gmm::sub_vector(rrhs, I2));
2612 if (brick.pbr->is_linear())
2613 brick.terms_to_be_computed =
false;
2625 if (version & BUILD_RHS) approx_external_load_ += brick.external_load;
2628 if (version & BUILD_RHS && version & BUILD_WITH_INTERNAL) {
2631 gmm::sub_interval IP(0,gmm::vect_size(rrhs));
2632 gmm::fill(full_rrhs, 0.);
2633 gmm::copy(rrhs, gmm::sub_vector(full_rrhs, IP));
2637 if (generic_expressions.size()) {
2640 if (version & BUILD_RHS)
2641 GMM_TRACE2(
"Global generic assembly RHS");
2642 if (version & BUILD_MATRIX)
2643 GMM_TRACE2(
"Global generic assembly tangent term");
2646 auto add_assignments_and_expressions_to_workspace =
2647 [&](ga_workspace &workspace)
2649 for (
const auto &ad : assignments)
2650 workspace.add_assignment_expression
2651 (ad.varname, ad.expr, ad.region, ad.order, ad.before);
2652 for (
const auto &ge : generic_expressions)
2653 workspace.add_expression(ge.expr, ge.mim, ge.region,
2654 2, ge.secondary_domain);
2657 const bool with_internal = version & BUILD_WITH_INTERNAL
2659 model_real_sparse_matrix intern_mat;
2660 model_real_plain_vector res0,
2663 size_type full_size = gmm::vect_size(full_rrhs),
2664 primary_size = gmm::vect_size(rrhs);
2666 if ((version & BUILD_RHS) || (version & BUILD_MATRIX && with_internal))
2667 gmm::resize(res0, with_internal ? full_size : primary_size);
2668 if (version & BUILD_MATRIX && with_internal)
2669 gmm::resize(res1, full_size);
2671 if (version & BUILD_MATRIX) {
2672 if (with_internal) {
2673 gmm::resize(intern_mat, full_size, primary_size);
2674 gmm::resize(res1, full_size);
2680 if (version & BUILD_RHS) {
2683 ga_workspace workspace(*
this);
2684 add_assignments_and_expressions_to_workspace(workspace);
2685 workspace.set_assembled_vector(res0_distro);
2686 workspace.assembly(1, with_internal);
2687 if (with_internal) {
2688 gmm::copy(res0_distro.get(), res1_distro.get());
2689 gmm::add(gmm::scaled(full_rrhs, scalar_type(-1)),
2691 workspace.set_assembled_vector(res1_distro);
2692 workspace.set_internal_coupling_matrix(intern_mat_distro);
2694 workspace.set_assembled_matrix(tangent_matrix_distro);
2695 workspace.assembly(2, with_internal);
2700 ga_workspace workspace(*
this);
2701 add_assignments_and_expressions_to_workspace(workspace);
2702 if (with_internal) {
2703 gmm::copy(gmm::scaled(full_rrhs, scalar_type(-1)),
2705 workspace.set_assembled_vector(res1_distro);
2706 workspace.set_internal_coupling_matrix(intern_mat_distro);
2708 workspace.set_assembled_matrix(tangent_matrix_distro);
2709 workspace.assembly(2, with_internal);
2713 else if (version & BUILD_RHS) {
2716 ga_workspace workspace(*
this);
2717 add_assignments_and_expressions_to_workspace(workspace);
2718 workspace.set_assembled_vector(res0_distro);
2719 workspace.assembly(1, with_internal);
2723 if (version & BUILD_RHS) {
2724 gmm::scale(res0, scalar_type(-1));
2725 if (with_internal) {
2726 gmm::sub_interval IP(0,gmm::vect_size(rrhs));
2727 gmm::add(gmm::sub_vector(res0, IP), rrhs);
2728 gmm::add(res0, full_rrhs);
2730 gmm::add(res0, rrhs);
2733 if (version & BUILD_MATRIX && with_internal) {
2734 gmm::scale(res1, scalar_type(-1));
2735 gmm::sub_interval IP(0, primary_size),
2736 II(primary_size, full_size-primary_size);
2737 gmm::copy(gmm::sub_matrix(intern_mat, II, IP), internal_rTM);
2738 gmm::add(gmm::sub_vector(res1, IP), rrhs);
2739 gmm::copy(gmm::sub_vector(res1, II), internal_sol);
2744 if ((version & BUILD_RHS) || (version & BUILD_MATRIX)) {
2746 std::vector<size_type> dof_indices;
2747 std::vector<complex_type> dof_pr_values;
2748 std::vector<complex_type> dof_go_values;
2749 for (
const auto &keyval : complex_dof_constraints) {
2750 const gmm::sub_interval &I = interval_of_variable(keyval.first);
2752 for (
const auto &val : keyval.second) {
2753 dof_indices.push_back(val.first + I.first());
2754 dof_go_values.push_back(val.second);
2755 dof_pr_values.push_back(V[val.first]);
2759 if (dof_indices.size()) {
2760 gmm::sub_index SI(dof_indices);
2761 gmm::sub_interval II(0,
nb_dof());
2763 if (version & BUILD_RHS) {
2764 if (MPI_IS_MASTER())
2765 approx_external_load_ += gmm::vect_norm1(dof_go_values);
2767 if (is_symmetric_) {
2768 scalar_type valnorm = gmm::vect_norm2(dof_go_values);
2769 if (valnorm > scalar_type(0)) {
2770 GMM_ASSERT1(version & BUILD_MATRIX,
"Rhs only for a "
2771 "symmetric linear problem with dof "
2772 "constraint not allowed");
2773 model_complex_plain_vector vv(gmm::vect_size(crhs));
2774 gmm::mult(gmm::sub_matrix(cTM, II, SI), dof_go_values, vv);
2776 gmm::add(gmm::scaled(vv, scalar_type(-1)), crhs);
2779 gmm::copy(dof_go_values, gmm::sub_vector(crhs, SI));
2781 gmm::add(dof_go_values,
2782 gmm::scaled(dof_pr_values, complex_type(-1)),
2783 gmm::sub_vector(crhs, SI));
2786 if (version & BUILD_MATRIX) {
2787 gmm::clear(gmm::sub_matrix(cTM, SI, II));
2788 if (is_symmetric_) gmm::clear(gmm::sub_matrix(cTM, II, SI));
2790 if (MPI_IS_MASTER()) {
2791 for (
size_type i = 0; i < dof_indices.size(); ++i)
2792 cTM(dof_indices[i], dof_indices[i]) = complex_type(1);
2797 std::vector<size_type> dof_indices;
2798 std::vector<scalar_type> dof_pr_values;
2799 std::vector<scalar_type> dof_go_values;
2800 for (
const auto &keyval : real_dof_constraints) {
2801 const gmm::sub_interval &I = interval_of_variable(keyval.first);
2802 const model_real_plain_vector &V =
real_variable(keyval.first);
2803 for (
const auto &val : keyval.second) {
2804 dof_indices.push_back(val.first + I.first());
2805 dof_go_values.push_back(val.second);
2806 dof_pr_values.push_back(V[val.first]);
2810 #if GETFEM_PARA_LEVEL > 1
2811 GMM_ASSERT1(MPI_IS_MASTER() || (dof_indices.size() == 0),
2812 "Sorry, for the moment, the dof constraints have to be "
2813 "added on the master process only");
2814 size_type dof_indices_size = dof_indices.size();
2815 MPI_BCAST0_SCALAR(dof_indices_size);
2816 dof_indices.resize(dof_indices_size);
2817 MPI_BCAST0_VECTOR(dof_indices);
2818 dof_pr_values.resize(dof_indices_size);
2819 MPI_BCAST0_VECTOR(dof_pr_values);
2820 dof_go_values.resize(dof_indices_size);
2821 MPI_BCAST0_VECTOR(dof_go_values);
2824 if (dof_indices.size()) {
2825 gmm::sub_index SI(dof_indices);
2826 gmm::sub_interval II(0,
nb_dof());
2828 if (version & BUILD_RHS) {
2829 if (MPI_IS_MASTER())
2830 approx_external_load_ += gmm::vect_norm1(dof_go_values);
2832 if (is_symmetric_) {
2833 scalar_type valnorm = gmm::vect_norm2(dof_go_values);
2834 if (valnorm > scalar_type(0)) {
2835 GMM_ASSERT1(version & BUILD_MATRIX,
"Rhs only for a "
2836 "symmetric linear problem with dof "
2837 "constraint not allowed");
2838 model_real_plain_vector vv(gmm::vect_size(rrhs));
2839 gmm::mult(gmm::sub_matrix(rTM, II, SI), dof_go_values, vv);
2841 gmm::add(gmm::scaled(vv, scalar_type(-1)), rrhs);
2844 gmm::copy(dof_go_values, gmm::sub_vector(rrhs, SI));
2846 gmm::add(dof_go_values,
2847 gmm::scaled(dof_pr_values, scalar_type(-1)),
2848 gmm::sub_vector(rrhs, SI));
2851 if (version & BUILD_MATRIX) {
2852 gmm::clear(gmm::sub_matrix(rTM, SI, II));
2853 if (is_symmetric_) gmm::clear(gmm::sub_matrix(rTM, II, SI));
2855 if (MPI_IS_MASTER()) {
2856 for (
size_type i = 0; i < dof_indices.size(); ++i)
2857 rTM(dof_indices[i], dof_indices[i]) = scalar_type(1);
2864 if (version & BUILD_RHS) {
2865 approx_external_load_ = MPI_SUM_SCALAR(approx_external_load_);
2868 #if GETFEM_PARA_LEVEL > 1
2870 if (MPI_IS_MASTER()) cout <<
"Assembly time " << MPI_Wtime()-t_ref << endl;
2879 return it->second.associated_mf();
2885 return it->second.passociated_mf();
2889 model::qdims_of_variable(
const std::string &name)
const {
2891 const mesh_fem *mf = it->second.passociated_mf();
2892 const im_data *imd = it->second.imd;
2895 bgeot::multi_index mi = mf->get_qdims();
2896 if (n > 1 || it->second.qdims.size() > 1) {
2898 if (mi.back() == 1) { mi.back() *= it->second.qdims[0]; ++i; }
2899 for (; i < it->second.qdims.size(); ++i)
2900 mi.push_back(it->second.qdims[i]);
2904 bgeot::multi_index mi = imd->tensor_size();
2907 "Invalid mesh im data vector");
2908 if (n > 1 || it->second.qdims.size() > 1) {
2910 if (mi.back() == 1) { mi.back() *= it->second.qdims[0]; ++i; }
2911 for (; i < it->second.qdims.size(); ++i)
2912 mi.push_back(it->second.qdims[i]);
2916 return it->second.qdims;
2919 size_type model::qdim_of_variable(
const std::string &name)
const {
2921 const mesh_fem *mf = it->second.passociated_mf();
2922 const im_data *imd = it->second.imd;
2925 return mf->get_qdim() * n;
2927 return imd->tensor_size().total_size() * n;
2933 const gmm::sub_interval &
2934 model::interval_of_variable(
const std::string &name)
const {
2936 if (act_size_to_be_done) actualize_sizes();
2937 VAR_SET::const_iterator it = find_variable(name);
2938 return it->second.I;
2941 const model_real_plain_vector &
2947 const model_real_plain_vector &
2949 GMM_ASSERT1(!complex_version,
"This model is a complex one");
2950 GMM_ASSERT1(!
is_old(name),
"Please don't use Old_ prefix in combination "
2951 "with variable version");
2953 auto it = variables.find(name);
2954 GMM_ASSERT1(it != variables.end(),
"Undefined variable " << name);
2955 if (act_size_to_be_done && it->second.is_fem_dofs) {
2956 if (it->second.filter != VDESCRFILTER_NO)
2959 it->second.set_size();
2961 if (niter ==
size_type(-1)) niter = it->second.default_iter;
2962 GMM_ASSERT1(it->second.n_iter + it->second.n_temp_iter > niter,
2963 "Invalid iteration number " << niter <<
" for " << name);
2964 return it->second.real_value[niter];
2967 const model_complex_plain_vector &
2973 const model_complex_plain_vector &
2975 GMM_ASSERT1(complex_version,
"This model is a real one");
2976 GMM_ASSERT1(!
is_old(name),
"Please don't use Old_ prefix in combination with"
2977 " variable version");
2979 auto it = variables.find(name);
2980 GMM_ASSERT1(it!=variables.end(),
"Undefined variable " << name);
2981 if (act_size_to_be_done && it->second.is_fem_dofs) {
2982 if (it->second.filter != VDESCRFILTER_NO)
2985 it->second.set_size();
2987 if (niter ==
size_type(-1)) niter = it->second.default_iter;
2988 GMM_ASSERT1(it->second.n_iter + it->second.n_temp_iter > niter,
2989 "Invalid iteration number "
2990 << niter <<
" for " << name);
2991 return it->second.complex_value[niter];
2994 model_real_plain_vector &
3001 model_real_plain_vector &
3003 GMM_ASSERT1(!complex_version,
"This model is a complex one");
3004 GMM_ASSERT1(!
is_old(name),
"Please don't use Old_ prefix in combination with"
3005 " variable version");
3007 auto it = variables.find(name);
3008 GMM_ASSERT1(it!=variables.end(),
"Undefined variable " << name);
3009 if (act_size_to_be_done && it->second.is_fem_dofs) {
3010 if (it->second.filter != VDESCRFILTER_NO)
3013 it->second.set_size();
3015 if (niter ==
size_type(-1)) niter = it->second.default_iter;
3016 it->second.v_num_data[niter] = act_counter();
3017 GMM_ASSERT1(it->second.n_iter + it->second.n_temp_iter > niter,
3018 "Invalid iteration number "
3019 << niter <<
" for " << name);
3020 return it->second.real_value[niter];
3023 model_complex_plain_vector &
3029 model_complex_plain_vector &
3031 GMM_ASSERT1(complex_version,
"This model is a real one");
3032 GMM_ASSERT1(!
is_old(name),
"Please don't use Old_ prefix in combination with"
3033 " variable version");
3035 auto it = variables.find(name);
3036 GMM_ASSERT1(it!=variables.end(),
"Undefined variable " << name);
3037 if (act_size_to_be_done && it->second.is_fem_dofs) {
3038 if (it->second.filter != VDESCRFILTER_NO)
3041 it->second.set_size();
3043 if (niter ==
size_type(-1)) niter = it->second.default_iter;
3044 it->second.v_num_data[niter] = act_counter();
3045 GMM_ASSERT1(it->second.n_iter + it->second.n_temp_iter > niter,
3046 "Invalid iteration number "
3047 << niter <<
" for " << name);
3048 return it->second.complex_value[niter];
3051 model_real_plain_vector &
3052 model::set_real_constant_part(
const std::string &name)
const {
3053 GMM_ASSERT1(!complex_version,
"This model is a complex one");
3055 VAR_SET::iterator it = variables.find(name);
3056 GMM_ASSERT1(it != variables.end(),
"Undefined variable " << name);
3057 GMM_ASSERT1(it->second.is_affine_dependent,
3058 "Only for affine dependent variables");
3059 if (act_size_to_be_done && it->second.is_fem_dofs) {
3060 if (it->second.filter != VDESCRFILTER_NO)
3063 it->second.set_size();
3065 for (
auto &v_num : it->second.v_num_data) v_num = act_counter();
3066 return it->second.affine_real_value;
3069 model_complex_plain_vector &
3070 model::set_complex_constant_part(
const std::string &name)
const {
3071 GMM_ASSERT1(complex_version,
"This model is a real one");
3073 VAR_SET::iterator it = variables.find(name);
3074 GMM_ASSERT1(it!=variables.end(),
"Undefined variable " << name);
3075 if (act_size_to_be_done && it->second.is_fem_dofs) {
3076 if (it->second.filter != VDESCRFILTER_NO)
3079 it->second.set_size();
3081 for (
auto &v_num : it->second.v_num_data) v_num = act_counter();
3082 return it->second.affine_complex_value;
3087 const brick_description &brick = bricks[ind_brick];
3088 update_brick(ind_brick, model::BUILD_ALL);
3090 brick.pbr->check_stiffness_matrix_and_rhs(*
this, ind_brick, brick.tlist,
3091 brick.vlist, brick.dlist, brick.mims, brick.rmatlist,
3092 brick.rveclist[0], brick.rveclist_sym[0], brick.region);
3095 void model::clear() {
3097 active_bricks.clear();
3098 valid_bricks.clear();
3099 real_dof_constraints.clear();
3100 complex_dof_constraints.clear();
3102 rTM = model_real_sparse_matrix();
3103 cTM = model_complex_sparse_matrix();
3104 rrhs = model_real_plain_vector();
3105 crhs = model_complex_plain_vector();
3110 void virtual_brick::full_asm_real_tangent_terms_(
const model &md,
size_type ind_brick,
3111 const model::varnamelist &term_list,
3112 const model::varnamelist &var_list,
3113 const model::mimlist &mim_list,
3114 model::real_matlist &rmatlist,
3115 model::real_veclist &rveclist,
3116 model::real_veclist &rveclist_sym,
3117 size_type region, build_version version)
const
3120 rveclist, rveclist_sym, region, version);
3122 rveclist, rveclist_sym, region, version);
3124 rveclist, rveclist_sym, region, version);
3129 const model::termlist& tlist,
3130 const model::varnamelist &vl,
3131 const model::varnamelist &dl,
3132 const model::mimlist &mims,
3133 model::real_matlist &matl,
3134 model::real_veclist &rvc1,
3135 model::real_veclist &rvc2,
3137 const scalar_type TINY)
const
3139 std::cout<<
"******Verifying stiffnesses of *******"<<std::endl;
3140 std::cout<<
"*** "<<brick_name()<<std::endl;
3143 std::map<std::string,size_type> rhs_index;
3144 for(
size_type iterm=0;iterm<matl.size();iterm++)
3145 if (tlist[iterm].var1==tlist[iterm].var2) rhs_index[tlist[iterm].var1]=iterm;
3147 if (rhs_index.size()==0){
3148 GMM_WARNING0(
"*** cannot verify stiffness for this brick***");
3151 full_asm_real_tangent_terms_(md, s, vl, dl, mims, matl, rvc1, rvc2,
3152 rg, model::BUILD_MATRIX);
3153 for(
size_type iterm=0;iterm<matl.size();iterm++)
3156 std::cout<<std::endl;
3157 std::cout<<
" Stiffness["<<tlist[iterm].var1
3158 <<
","<<tlist[iterm].var2<<
"]:"<<std::endl;
3161 std::cout<<
" "<<tlist[iterm].var1<<
" has zero size. Skipping this term"<<std::endl;
3166 std::cout<<
" "<<tlist[iterm].var2<<
" has zero size. Skipping this term"<<std::endl;
3170 model_real_sparse_matrix SM(matl[iterm]);
3171 gmm::fill(rvc1[rhs_index[tlist[iterm].var1]], 0.0);
3172 full_asm_real_tangent_terms_(md, s, vl, dl, mims, matl, rvc1, rvc2,
3173 rg, model::BUILD_RHS);
3174 if (gmm::mat_euclidean_norm(matl[iterm])<1e-12){
3175 std::cout<<
" The assembled matrix is nearly zero, skipping."<<std::endl;
3178 model_real_plain_vector RHS0(rvc1[rhs_index[tlist[iterm].var1]]);
3181 model_real_sparse_matrix fdSM(matl[iterm].nrows(), matl[iterm].ncols());
3183 model_real_plain_vector& RHS1 = rvc1[rhs_index[tlist[iterm].var1]];
3185 scalar_type relative_tiny = gmm::vect_norminf(RHS1)*TINY;
3186 if (relative_tiny < TINY) relative_tiny = TINY;
3188 for (
size_type j = 0; j < matl[iterm].ncols(); j++){
3189 U[j] += relative_tiny;
3190 gmm::fill(RHS1, 0.0);
3191 full_asm_real_tangent_terms_(md, s, vl, dl, mims, matl, rvc1, rvc2,
3192 rg, model::BUILD_RHS);
3193 for (
size_type i = 0; i<matl[iterm].nrows(); i++)
3194 fdSM(i, j) = (RHS0[i] - RHS1[i]) / relative_tiny;
3195 U[j] -= relative_tiny;
3198 model_real_sparse_matrix diffSM(matl[iterm].nrows(),matl[iterm].ncols());
3199 gmm::add(SM,gmm::scaled(fdSM,-1.0),diffSM);
3200 scalar_type norm_error_euc
3201 = gmm::mat_euclidean_norm(diffSM)/gmm::mat_euclidean_norm(SM)*100;
3202 scalar_type norm_error_1
3203 = gmm::mat_norm1(diffSM)/gmm::mat_norm1(SM)*100;
3204 scalar_type norm_error_max
3205 = gmm::mat_maxnorm(diffSM)/gmm::mat_maxnorm(SM)*100;
3208 scalar_type nsym_norm_error_euc=0.0;
3209 scalar_type nsym_norm_error_1=0.0;
3210 scalar_type nsym_norm_error_max=0.0;
3211 if (tlist[iterm].var1==tlist[iterm].var2){
3212 model_real_sparse_matrix diffSMtransposed(matl[iterm].nrows(),matl[iterm].ncols());
3213 gmm::add(gmm::transposed(fdSM),gmm::scaled(fdSM,-1.0),diffSMtransposed);
3215 = gmm::mat_euclidean_norm(diffSMtransposed)/gmm::mat_euclidean_norm(fdSM)*100;
3217 = gmm::mat_norm1(diffSMtransposed)/gmm::mat_norm1(fdSM)*100;
3219 = gmm::mat_maxnorm(diffSMtransposed)/gmm::mat_maxnorm(fdSM)*100;
3223 if(rvc1[0].size()<8){
3224 std::cout <<
"RHS Stiffness Matrix: \n";
3225 std::cout <<
"------------------------\n";
3226 for(
size_type i=0; i < rvc1[iterm].size(); ++i){
3228 for(
size_type j=0; j < rvc1[iterm].size(); ++j){
3229 std::cout << fdSM(i,j) <<
" ";
3233 std::cout <<
"Analytical Stiffness Matrix: \n";
3234 std::cout <<
"------------------------\n";
3235 for(
size_type i=0; i < rvc1[iterm].size(); ++i){
3237 for(
size_type j=0; j < rvc1[iterm].size(); ++j){
3238 std::cout << matl[iterm](i,j) <<
" ";
3242 std::cout <<
"Vector U: \n";
3243 std::cout <<
"------------------------\n";
3244 for(
size_type i=0; i < rvc1[iterm].size(); ++i){
3251 <<
"\n\nfinite diff test error_norm_eucl: " << norm_error_euc <<
"%\n"
3252 <<
"finite diff test error_norm1: " << norm_error_1 <<
"%\n"
3253 <<
"finite diff test error_max_norm: " << norm_error_max <<
"%\n\n\n";
3255 if (tlist[iterm].var1==tlist[iterm].var2){
3257 <<
"Nonsymmetrical test error_norm_eucl: "<< nsym_norm_error_euc<<
"%\n"
3258 <<
"Nonsymmetrical test error_norm1: " << nsym_norm_error_1 <<
"%\n"
3259 <<
"Nonsymmetrical test error_max_norm: " << nsym_norm_error_max <<
"%"
3279 struct gen_source_term_assembly_brick :
public virtual_brick {
3281 std::string expr, directvarname, directdataname;
3282 model::varnamelist vl_test1;
3283 std::string secondary_domain;
3286 const model::varnamelist &,
3287 const model::varnamelist &,
3288 const model::mimlist &mims,
3289 model::real_matlist &,
3290 model::real_veclist &vecl,
3291 model::real_veclist &,
3293 build_version)
const override {
3294 GMM_ASSERT1(vecl.size() == vl_test1.size()
3295 + ((directdataname.size() == 0) ? 0 : 1),
"Wrong number "
3296 "of terms for Generic source term assembly brick ");
3297 GMM_ASSERT1(mims.size() == 1,
"Generic source term assembly brick "
3298 "needs one and only one mesh_im");
3299 GMM_TRACE2(
"Generic source term assembly");
3301 gmm::clear(vecl[0]);
3305 ga_workspace workspace(md, ga_workspace::inherit::ALL);
3306 GMM_TRACE2(name <<
": generic source term assembly");
3307 workspace.add_expression(expr, *(mims[0]), region, 1, secondary_domain);
3308 workspace.assembly(1);
3309 const auto &V=workspace.assembled_vector();
3310 for (
size_type i = 0; i < vl_test1.size(); ++i) {
3311 const auto &I=workspace.interval_of_variable(vl_test1[i]);
3312 gmm::copy(gmm::sub_vector(V, I), vecl[i]);
3316 if (directvarname.size()) {
3321 void real_post_assembly_in_serial(
const model &md,
size_type ib,
3322 const model::varnamelist &,
3323 const model::varnamelist &,
3324 const model::mimlist &,
3325 model::real_matlist &,
3326 model::real_veclist &vecl,
3327 model::real_veclist &,
3329 build_version)
const override {
3330 scalar_type el = scalar_type(0);
3331 for (
size_type i=0; i < vecl.size(); ++i) el += gmm::vect_norm1(vecl[i]);
3332 md.add_external_load(ib, el);
3335 virtual std::string declare_volume_assembly_string
3336 (
const model &,
size_type,
const model::varnamelist &,
3337 const model::varnamelist &)
const {
3338 return std::string();
3341 gen_source_term_assembly_brick(
const std::string &expr_,
3342 std::string brickname,
3343 const model::varnamelist &vl_test1_,
3344 const std::string &directvarname_,
3345 const std::string &directdataname_,
3346 const std::string &secdom)
3347 : vl_test1(vl_test1_), secondary_domain(secdom) {
3348 if (brickname.size() == 0)
3349 brickname =
"Generic source term assembly brick";
3351 set_flags(brickname,
true ,
3355 directvarname = directvarname_; directdataname = directdataname_;
3360 static bool check_compatibility_vl_test(model &md,
3361 const model::varnamelist vl_test) {
3362 model::varnamelist org;
3363 for (
size_type i = 0; i < vl_test.size(); ++i) {
3364 if (md.is_affine_dependent_variable(vl_test[i]))
3365 org.push_back(md.org_variable(vl_test[i]));
3367 for (
size_type i = 0; i < vl_test.size(); ++i)
3368 for (
size_type j = 0; j < org.size(); ++j)
3369 if (vl_test[i].compare(org[j]) == 0)
return false;
3374 (model &md,
const mesh_im &mim,
const std::string &expr,
size_type region,
3375 const std::string &brickname, std::string directvarname,
3376 const std::string &directdataname,
bool return_if_nonlin,
3377 const std::string &secondary_domain) {
3379 ga_workspace workspace(md);
3380 size_type order = workspace.add_expression(expr, mim, region, 1,
3382 GMM_ASSERT1(order <= 1,
"Wrong order for a source term");
3383 model::varnamelist vl, vl_test1, vl_test2, dl;
3384 bool is_lin = workspace.used_variables(vl, vl_test1, vl_test2, dl, 1);
3385 if (!is_lin && return_if_nonlin)
return size_type(-1);
3386 GMM_ASSERT1(is_lin,
"Nonlinear term");
3387 GMM_ASSERT1(check_compatibility_vl_test(md, vl_test1),
3388 "This brick do not support the assembly on both an affine "
3389 "dependent variable and its original variable. "
3390 "Split the brick.");
3392 if (directdataname.size()) {
3393 vl.push_back(directvarname);
3394 dl.push_back(directdataname);
3395 }
else directvarname =
"";
3397 pbrick pbr = std::make_shared<gen_source_term_assembly_brick>
3398 (expr, brickname, vl_test1, directvarname, directdataname,
3402 for (
size_type i = 0; i < vl_test1.size(); ++i)
3403 tl.push_back(model::term_description(vl_test1[i]));
3404 if (directdataname.size())
3405 tl.push_back(model::term_description(directvarname));
3407 return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
3412 const std::string &brickname,
const std::string &directvarname,
3413 const std::string &directdataname,
bool return_if_nonlin) {
3414 return add_source_term_(md, mim, expr, region, brickname, directvarname,
3415 directdataname, return_if_nonlin,
"");
3420 const std::string &secondary_domain,
3421 const std::string &brickname,
const std::string &directvarname,
3422 const std::string &directdataname,
bool return_if_nonlin) {
3423 return add_source_term_(md, mim, expr, region, brickname, directvarname,
3424 directdataname, return_if_nonlin, secondary_domain);
3433 struct gen_linear_assembly_brick :
public virtual_brick {
3437 model::varnamelist vl_test1, vl_test2;
3438 std::string secondary_domain;
3440 virtual void asm_real_tangent_terms(
const model &md,
size_type ib,
3441 const model::varnamelist &,
3442 const model::varnamelist &dl,
3443 const model::mimlist &mims,
3444 model::real_matlist &matl,
3445 model::real_veclist &,
3446 model::real_veclist &,
3448 build_version version)
const {
3449 GMM_ASSERT1(matl.size() == vl_test1.size(),
3450 "Wrong number of terms for Generic linear assembly brick");
3451 GMM_ASSERT1(mims.size() == 1,
3452 "Generic linear assembly brick needs one and only one "
3454 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0);
3455 for (
size_type i = 0; i < dl.size(); ++i)
3456 recompute_matrix = recompute_matrix ||
3457 md.is_var_newer_than_brick(dl[i], ib);
3459 if (recompute_matrix) {
3461 ga_workspace workspace(md, ga_workspace::inherit::ALL);
3462 workspace.add_expression(expr, *(mims[0]), region, 2, secondary_domain);
3463 GMM_TRACE2(name <<
": generic matrix assembly");
3464 workspace.assembly(2);
3465 const auto &R=workspace.assembled_matrix();
3466 for (
size_type i = 0; i < vl_test1.size(); ++i) {
3467 scalar_type alpha = scalar_type(1)
3468 / ( workspace.factor_of_variable(vl_test1[i]) *
3469 workspace.factor_of_variable(vl_test2[i]));
3470 const auto &I1=workspace.interval_of_variable(vl_test1[i]),
3471 &I2=workspace.interval_of_variable(vl_test2[i]);
3472 gmm::copy(gmm::scaled(gmm::sub_matrix(R, I1, I2), alpha),
3478 virtual std::string declare_volume_assembly_string
3479 (
const model &,
size_type,
const model::varnamelist &,
3480 const model::varnamelist &)
const {
3481 return is_lower_dim ? std::string() : expr;
3484 gen_linear_assembly_brick(
const std::string &expr_,
const mesh_im &mim,
3486 bool is_coer, std::string brickname,
3487 const model::varnamelist &vl_test1_,
3488 const model::varnamelist &vl_test2_,
3489 const std::string &secdom)
3490 : vl_test1(vl_test1_), vl_test2(vl_test2_), secondary_domain(secdom) {
3491 if (brickname.size() == 0) brickname =
"Generic linear assembly brick";
3493 is_lower_dim = mim.is_lower_dimensional();
3494 set_flags(brickname,
true ,
3501 static bool check_compatibility_vl_test(model &md,
3502 const model::varnamelist vl_test1,
3503 const model::varnamelist vl_test2) {
3504 for (
size_type i = 0; i < vl_test1.size(); ++i)
3505 for (
size_type j = 0; j < vl_test1.size(); ++j) {
3506 bool is1 = md.is_affine_dependent_variable(vl_test1[i]);
3507 bool is2 = md.is_affine_dependent_variable(vl_test2[i]);
3509 const std::string &org1
3510 = is1 ? md.org_variable(vl_test1[i]) : vl_test1[i];
3511 const std::string &org2
3512 = is2 ? md.org_variable(vl_test2[i]) : vl_test2[i];
3513 bool is1_bis = md.is_affine_dependent_variable(vl_test1[j]);
3514 bool is2_bis = md.is_affine_dependent_variable(vl_test2[j]);
3515 const std::string &org1_bis = is1_bis ? md.org_variable(vl_test1[j])
3517 const std::string &org2_bis = is2_bis ? md.org_variable(vl_test2[j])
3519 if (org1.compare(org1_bis) == 0 && org2.compare(org2_bis))
3529 (model &md,
const mesh_im &mim,
const std::string &expr,
size_type region,
3530 bool is_sym,
bool is_coercive,
const std::string &brickname,
3531 bool return_if_nonlin,
const std::string &secondary_domain) {
3533 ga_workspace workspace(md, ga_workspace::inherit::ALL);
3534 size_type order = workspace.add_expression(expr, mim, region,
3535 2, secondary_domain);
3536 model::varnamelist vl, vl_test1, vl_test2, dl;
3537 bool is_lin = workspace.used_variables(vl, vl_test1, vl_test2, dl, 2);
3539 if (!is_lin && return_if_nonlin)
return size_type(-1);
3540 GMM_ASSERT1(is_lin,
"Nonlinear term");
3541 if (order == 0) { is_coercive = is_sym =
true; }
3543 std::string const_expr= workspace.extract_constant_term(mim.linked_mesh());
3544 if (const_expr.size()) {
3546 (md, mim, const_expr, region, brickname+
" (source term)",
3547 "",
"",
false, secondary_domain);
3552 GMM_ASSERT1(check_compatibility_vl_test(md, vl_test1, vl_test2),
3553 "This brick do not support the assembly on both an affine "
3554 "dependent variable and its original variable. "
3555 "Split the brick.");
3557 if (vl_test1.size()) {
3558 pbrick pbr = std::make_shared<gen_linear_assembly_brick>
3559 (expr, mim, is_sym, is_coercive, brickname, vl_test1, vl_test2,
3562 for (
size_type i = 0; i < vl_test1.size(); ++i)
3563 tl.push_back(model::term_description(vl_test1[i], vl_test2[i],
false));
3565 return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
3572 bool is_sym,
bool is_coercive,
const std::string &brickname,
3573 bool return_if_nonlin) {
3574 return add_linear_term_(md, mim, expr, region, is_sym, is_coercive,
3575 brickname, return_if_nonlin,
"");
3580 const std::string &secondary_domain,
bool is_sym,
bool is_coercive,
3581 const std::string &brickname,
bool return_if_nonlin) {
3582 return add_linear_term_(md, mim, expr, region, is_sym, is_coercive,
3583 brickname, return_if_nonlin, secondary_domain);
3593 struct gen_nonlinear_assembly_brick :
public virtual_brick {
3597 std::string secondary_domain;
3599 virtual void real_post_assembly_in_serial(
const model &md,
size_type ,
3600 const model::varnamelist &,
3601 const model::varnamelist &,
3602 const model::mimlist &mims,
3603 model::real_matlist &,
3604 model::real_veclist &,
3605 model::real_veclist &,
3607 build_version)
const {
3608 GMM_ASSERT1(mims.size() == 1,
3609 "Generic linear assembly brick needs one and only one "
3611 md.add_generic_expression(expr, *(mims[0]), region, secondary_domain);
3614 virtual std::string declare_volume_assembly_string
3615 (
const model &,
size_type,
const model::varnamelist &,
3616 const model::varnamelist &)
const {
3621 gen_nonlinear_assembly_brick(
const std::string &expr_,
const mesh_im &mim,
3624 std::string brickname,
3625 const std::string &secdom) {
3626 if (brickname.size() == 0) brickname =
"Generic linear assembly brick";
3628 secondary_domain = secdom;
3629 is_lower_dim = mim.is_lower_dimensional();
3630 set_flags(brickname,
false ,
3638 (model &md,
const mesh_im &mim,
const std::string &expr,
size_type region,
3639 bool is_sym,
bool is_coercive,
const std::string &brickname,
3640 const std::string &secondary_domain) {
3642 ga_workspace workspace(md);
3643 size_type order = workspace.add_expression(expr, mim, region, 2,
3645 GMM_ASSERT1(order < 2,
"Order two test functions (Test2) are not allowed"
3646 " in assembly string for nonlinear terms");
3647 model::varnamelist vl, vl_test1, vl_test2, ddl, dl;
3648 workspace.used_variables(vl, vl_test1, vl_test2, ddl, order);
3649 for (
size_type i = 0; i < ddl.size(); ++i)
3650 if (md.is_true_data(ddl[i])) dl.push_back(ddl[i]);
3651 else vl.push_back(ddl[i]);
3652 if (order == 0) { is_coercive = is_sym =
true; }
3653 pbrick pbr = std::make_shared<gen_nonlinear_assembly_brick>
3654 (expr, mim, is_sym, is_coercive, brickname, secondary_domain);
3658 return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
3664 bool is_sym,
bool is_coercive,
const std::string &brickname) {
3665 return add_nonlinear_term_(md, mim, expr, region, is_sym, is_coercive,
3671 const std::string &secondary_domain,
bool is_sym,
bool is_coercive,
3672 const std::string &brickname) {
3673 return add_nonlinear_term_(md, mim, expr, region, is_sym, is_coercive,
3674 brickname, secondary_domain);
3685 struct generic_elliptic_brick :
public virtual_brick {
3687 virtual void asm_real_tangent_terms(
const model &md,
size_type ,
3688 const model::varnamelist &vl,
3689 const model::varnamelist &dl,
3690 const model::mimlist &mims,
3691 model::real_matlist &matl,
3692 model::real_veclist &,
3693 model::real_veclist &,
3695 build_version)
const {
3696 GMM_ASSERT1(matl.size() == 1,
3697 "Generic elliptic brick has one and only one term");
3698 GMM_ASSERT1(mims.size() == 1,
3699 "Generic elliptic brick need one and only one mesh_im");
3700 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
3701 "Wrong number of variables for generic elliptic brick");
3703 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
3704 const mesh &m = mf_u.linked_mesh();
3705 size_type N = m.dim(), Q = mf_u.get_qdim(), s = 1;
3706 const mesh_im &mim = *mims[0];
3707 const model_real_plain_vector *A = 0;
3708 const mesh_fem *mf_a = 0;
3709 mesh_region rg(region);
3710 m.intersect_with_mpi_region(rg);
3711 if (dl.size() > 0) {
3712 A = &(md.real_variable(dl[0]));
3713 mf_a = md.pmesh_fem_of_variable(dl[0]);
3714 s = gmm::vect_size(*A);
3715 if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
3719 GMM_TRACE2(
"Generic elliptic term assembly");
3724 (matl[0], mim, mf_u, *mf_a, *A, rg);
3727 (matl[0], mim, mf_u, *mf_a, *A, rg);
3732 (matl[0], mim, mf_u, rg);
3735 (matl[0], mim, mf_u, rg);
3736 if (A) gmm::scale(matl[0], (*A)[0]);
3738 }
else if (s == N*N) {
3742 (matl[0], mim, mf_u, *mf_a, *A, rg);
3745 (matl[0], mim, mf_u, *mf_a, *A, rg);
3749 (matl[0], mim, mf_u, *A, rg);
3752 (matl[0], mim, mf_u, *A, rg);
3754 }
else if (s == N*N*Q*Q) {
3757 (matl[0], mim, mf_u, *mf_a, *A, rg);
3760 (matl[0], mim, mf_u, *A, rg);
3762 GMM_ASSERT1(
false,
"Bad format generic elliptic brick coefficient");
3766 virtual void real_post_assembly_in_serial(
const model &md,
size_type,
3767 const model::varnamelist &,
3768 const model::varnamelist &dl,
3769 const model::mimlist &,
3770 model::real_matlist &,
3771 model::real_veclist &,
3772 model::real_veclist &,
3774 build_version)
const
3776 const model_real_plain_vector *A = 0;
3777 const mesh_fem *mf_a = 0;
3780 A = &(md.real_variable(dl[0]));
3781 mf_a = md.pmesh_fem_of_variable(dl[0]);
3785 virtual void complex_post_assembly_in_serial(
const model &md,
size_type,
3786 const model::varnamelist &,
3787 const model::varnamelist &dl,
3788 const model::mimlist &,
3789 model::complex_matlist &,
3790 model::complex_veclist &,
3791 model::complex_veclist &,
3793 build_version)
const
3795 const model_real_plain_vector *A = 0;
3796 const mesh_fem *mf_a = 0;
3799 A = &(md.real_variable(dl[0]));
3800 mf_a = md.pmesh_fem_of_variable(dl[0]);
3804 virtual scalar_type asm_complex_tangent_terms(
const model &md,
size_type,
3805 const model::varnamelist &vl,
3806 const model::varnamelist &,
3807 const model::mimlist &,
3808 model::complex_matlist &matl,
3809 model::complex_veclist &,
3810 model::complex_veclist &,
3812 const model_complex_plain_vector &U = md.complex_variable(vl[0]);
3813 return gmm::abs(gmm::vect_hp(matl[0], U, U)) / scalar_type(2);
3817 virtual void asm_complex_tangent_terms(
const model &md,
size_type,
3818 const model::varnamelist &vl,
3819 const model::varnamelist &dl,
3820 const model::mimlist &mims,
3821 model::complex_matlist &matl,
3822 model::complex_veclist &,
3823 model::complex_veclist &,
3825 build_version)
const {
3826 GMM_ASSERT1(matl.size() == 1,
3827 "Generic elliptic brick has one and only one term");
3828 GMM_ASSERT1(mims.size() == 1,
3829 "Generic elliptic brick need one and only one mesh_im");
3830 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
3831 "Wrong number of variables for generic elliptic brick");
3833 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
3834 const mesh &m = mf_u.linked_mesh();
3835 size_type N = m.dim(), Q = mf_u.get_qdim(), s = 1;
3836 const mesh_im &mim = *mims[0];
3837 const model_real_plain_vector *A = 0;
3838 const mesh_fem *mf_a = 0;
3839 mesh_region rg(region);
3840 m.intersect_with_mpi_region(rg);
3843 if (dl.size() > 0) {
3844 A = &(md.real_variable(dl[0]));
3845 mf_a = md.pmesh_fem_of_variable(dl[0]);
3846 s = gmm::vect_size(*A);
3847 if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
3851 GMM_TRACE2(
"Generic elliptic term assembly");
3856 (matl[0], mim, mf_u, *mf_a, *A, rg);
3859 (matl[0], mim, mf_u, *mf_a, *A, rg);
3864 (gmm::real_part(matl[0]), mim, mf_u, rg);
3867 (gmm::real_part(matl[0]), mim, mf_u, rg);
3868 if (A) gmm::scale(matl[0], (*A)[0]);
3870 }
else if (s == N*N) {
3874 (matl[0], mim, mf_u, *mf_a, *A, rg);
3877 (matl[0], mim, mf_u, *mf_a, *A, rg);
3881 (matl[0], mim, mf_u, *A, rg);
3884 (matl[0], mim, mf_u, *A, rg);
3886 }
else if (s == N*N*Q*Q) {
3889 (matl[0], mim, mf_u, *mf_a, *A, rg);
3892 (matl[0], mim, mf_u, *A, rg);
3895 "Bad format generic elliptic brick coefficient");
3898 generic_elliptic_brick() {
3899 set_flags(
"Generic elliptic",
true ,
3907 const std::string &varname,
3910 pbrick pbr = std::make_shared<generic_elliptic_brick>();
3912 tl.push_back(model::term_description(varname, varname,
true));
3913 return md.
add_brick(pbr, model::varnamelist(1, varname),
3914 model::varnamelist(), tl, model::mimlist(1, &mim),
3917 std::string test_varname
3918 =
"Test_" + sup_previous_and_dot_to_varname(varname);
3923 expr =
"Grad_"+varname+
".Grad_"+test_varname;
3925 expr =
"Grad_"+varname+
":Grad_"+test_varname;
3927 "Laplacian",
false);
3932 const std::string &varname,
3933 const std::string &dataname,
3936 pbrick pbr = std::make_shared<generic_elliptic_brick>();
3938 tl.push_back(model::term_description(varname, varname,
true));
3939 return md.
add_brick(pbr, model::varnamelist(1, varname),
3940 model::varnamelist(1, dataname), tl,
3941 model::mimlist(1, &mim), region);
3943 std::string test_varname
3944 =
"Test_" + sup_previous_and_dot_to_varname(varname);
3957 if (qdim_data != 1) {
3958 GMM_ASSERT1(qdim_data == gmm::sqr(dim),
3959 "Wrong data size for generic elliptic brick");
3960 expr =
"((Reshape("+dataname+
",meshdim,meshdim))*Grad_"+varname+
").Grad_"
3963 expr =
"(("+dataname+
")*Grad_"+varname+
").Grad_"+test_varname;
3966 if (qdim_data != 1) {
3967 if (qdim_data == gmm::sqr(dim))
3968 expr =
"((Reshape("+dataname+
",meshdim,meshdim))*Grad_"+varname+
"):Grad_"
3970 else if (qdim_data == gmm::sqr(gmm::sqr(dim)))
3971 expr =
"((Reshape("+dataname+
",meshdim,meshdim,meshdim,meshdim))*Grad_"
3972 +varname+
"):Grad_"+test_varname;
3973 else GMM_ASSERT1(
false,
"Wrong data size for generic elliptic brick");
3975 expr =
"(("+dataname+
")*Grad_"+varname+
"):Grad_"+test_varname;
3979 (md, mim, expr, region,
true,
true,
"Generic elliptic",
true);
3982 "Generic elliptic (nonlinear)");
3994 struct source_term_brick :
public virtual_brick {
3996 void asm_real_tangent_terms(
const model &md,
size_type ,
3997 const model::varnamelist &vl,
3998 const model::varnamelist &dl,
3999 const model::mimlist &mims,
4000 model::real_matlist &,
4001 model::real_veclist &vecl,
4002 model::real_veclist &,
4004 build_version)
const override {
4005 GMM_ASSERT1(vecl.size() == 1,
4006 "Source term brick has one and only one term");
4007 GMM_ASSERT1(mims.size() == 1,
4008 "Source term brick need one and only one mesh_im");
4009 GMM_ASSERT1(vl.size() == 1 && dl.size() > 0 && dl.size() <= 2,
4010 "Wrong number of variables for source term brick");
4012 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4013 const mesh_im &mim = *mims[0];
4014 const model_real_plain_vector &A = md.real_variable(dl[0]);
4015 const mesh_fem *mf_data = md.pmesh_fem_of_variable(dl[0]);
4016 mesh_region rg(region);
4017 mim.linked_mesh().intersect_with_mpi_region(rg);
4020 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4022 GMM_ASSERT1(mf_u.get_qdim() == s,
4023 dl[0] <<
": bad format of source term data. "
4024 "Detected dimension is " << s <<
" should be "
4027 GMM_TRACE2(
"Source term assembly");
4033 if (dl.size() > 1) gmm::add(md.real_variable(dl[1]), vecl[0]);
4036 void real_post_assembly_in_serial(
const model &md,
size_type ib,
4037 const model::varnamelist &,
4038 const model::varnamelist &,
4039 const model::mimlist &,
4040 model::real_matlist &,
4041 model::real_veclist &vecl,
4042 model::real_veclist &,
4044 build_version)
const override
4046 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4050 void asm_complex_tangent_terms(
const model &md,
size_type ,
4051 const model::varnamelist &vl,
4052 const model::varnamelist &dl,
4053 const model::mimlist &mims,
4054 model::complex_matlist &,
4055 model::complex_veclist &vecl,
4056 model::complex_veclist &,
4058 build_version)
const override {
4059 GMM_ASSERT1(vecl.size() == 1,
4060 "Source term brick has one and only one term");
4061 GMM_ASSERT1(mims.size() == 1,
4062 "Source term brick need one and only one mesh_im");
4063 GMM_ASSERT1(vl.size() == 1 && dl.size() > 0 && dl.size() <= 2,
4064 "Wrong number of variables for source term brick");
4066 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4067 const mesh_im &mim = *mims[0];
4068 const model_complex_plain_vector &A = md.complex_variable(dl[0]);
4069 const mesh_fem *mf_data = md.pmesh_fem_of_variable(dl[0]);
4070 mesh_region rg(region);
4071 mim.linked_mesh().intersect_with_mpi_region(rg);
4074 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4076 GMM_ASSERT1(mf_u.get_qdim() == s,
"Bad format of source term data");
4078 GMM_TRACE2(
"Source term assembly");
4084 if (dl.size() > 1) gmm::add(md.complex_variable(dl[1]), vecl[0]);
4087 void complex_post_assembly_in_serial(
const model &md,
4089 const model::varnamelist &,
4090 const model::varnamelist &,
4091 const model::mimlist &,
4092 model::complex_matlist &,
4093 model::complex_veclist &vecl,
4094 model::complex_veclist &,
4095 size_type, build_version)
const override
4097 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4102 source_term_brick() {
4103 set_flags(
"Source term",
true ,
4113 const std::string &varname,
4114 const std::string &dataexpr,
4116 const std::string &directdataname) {
4118 pbrick pbr = std::make_shared<source_term_brick>();
4120 tl.push_back(model::term_description(varname));
4121 model::varnamelist vdata(1, dataexpr);
4122 if (directdataname.size()) vdata.push_back(directdataname);
4123 return md.
add_brick(pbr, model::varnamelist(1, varname),
4124 vdata, tl, model::mimlist(1, &mim), region);
4126 std::string test_varname
4127 =
"Test_" + sup_previous_and_dot_to_varname(varname);
4132 expr =
"("+dataexpr+
")*"+test_varname;
4134 expr =
"("+dataexpr+
")."+test_varname;
4135 size_type ib = add_source_term_generic_assembly_brick
4136 (md, mim, expr, region,
"Source term", varname, directdataname,
true);
4139 "Source term (nonlinear)");
4140 if (directdataname.size())
4141 add_source_term_generic_assembly_brick
4142 (md, mim,
"", region,
"Source term", varname, directdataname);
4154 struct normal_source_term_brick :
public virtual_brick {
4156 void asm_real_tangent_terms(
const model &md,
size_type ,
4157 const model::varnamelist &vl,
4158 const model::varnamelist &dl,
4159 const model::mimlist &mims,
4160 model::real_matlist &,
4161 model::real_veclist &vecl,
4162 model::real_veclist &,
4164 build_version)
const override {
4165 GMM_ASSERT1(vecl.size() == 1,
4166 "Source term brick has one and only one term");
4167 GMM_ASSERT1(mims.size() == 1,
4168 "Source term brick need one and only one mesh_im");
4169 GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
4170 "Wrong number of variables for source term brick");
4172 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4173 const mesh_im &mim = *mims[0];
4174 const model_real_plain_vector &A = md.real_variable(dl[0]);
4175 const mesh_fem *mf_data = md.pmesh_fem_of_variable(dl[0]);
4176 mesh_region rg(region);
4177 mim.linked_mesh().intersect_with_mpi_region(rg);
4179 size_type s = gmm::vect_size(A), N = mf_u.linked_mesh().dim();
4180 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4182 GMM_ASSERT1(mf_u.get_qdim()*N == s,
4183 dl[0] <<
": bad format of normal source term data. "
4184 "Detected dimension is " << s <<
" should be "
4187 GMM_TRACE2(
"source term assembly");
4194 void real_post_assembly_in_serial(
const model &md,
size_type ib,
4195 const model::varnamelist &,
4196 const model::varnamelist &,
4197 const model::mimlist &,
4198 model::real_matlist &,
4199 model::real_veclist &vecl,
4200 model::real_veclist &,
4202 build_version)
const override {
4203 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4207 virtual void asm_complex_tangent_terms(
const model &md,
size_type ,
4208 const model::varnamelist &vl,
4209 const model::varnamelist &dl,
4210 const model::mimlist &mims,
4211 model::complex_matlist &,
4212 model::complex_veclist &vecl,
4213 model::complex_veclist &,
4215 build_version)
const {
4216 GMM_ASSERT1(vecl.size() == 1,
4217 "Source term brick has one and only one term");
4218 GMM_ASSERT1(mims.size() == 1,
4219 "Source term brick need one and only one mesh_im");
4220 GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
4221 "Wrong number of variables for source term brick");
4223 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4224 const mesh_im &mim = *mims[0];
4225 const model_complex_plain_vector &A = md.complex_variable(dl[0]);
4226 const mesh_fem *mf_data = md.pmesh_fem_of_variable(dl[0]);
4227 mesh_region rg(region);
4228 mim.linked_mesh().intersect_with_mpi_region(rg);
4230 size_type s = gmm::vect_size(A), N = mf_u.linked_mesh().dim();
4231 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4233 GMM_ASSERT1(s == mf_u.get_qdim()*N,
"Bad format of source term data");
4235 GMM_TRACE2(
"source term assembly");
4242 void complex_post_assembly_in_serial(
const model &md,
size_type ib,
4243 const model::varnamelist &,
4244 const model::varnamelist &,
4245 const model::mimlist &,
4246 model::complex_matlist &,
4247 model::complex_veclist &vecl,
4248 model::complex_veclist &,
4250 build_version)
const override {
4251 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4254 normal_source_term_brick() {
4255 set_flags(
"Normal source term",
true ,
4265 const std::string &varname,
4266 const std::string &dataexpr,
4269 pbrick pbr = std::make_shared<normal_source_term_brick>();
4271 tl.push_back(model::term_description(varname));
4272 model::varnamelist vdata(1, dataexpr);
4273 return md.
add_brick(pbr, model::varnamelist(1, varname),
4274 vdata, tl, model::mimlist(1, &mim), region);
4276 std::string test_varname
4277 =
"Test_" + sup_previous_and_dot_to_varname(varname);
4282 expr =
"(("+dataexpr+
").Normal)*"+test_varname;
4284 expr =
"(Reshape("+dataexpr+
",qdim("+varname
4285 +
"),meshdim)*Normal)."+test_varname;
4286 return add_source_term_generic_assembly_brick
4287 (md, mim, expr, region,
"Source term");
4300 struct Dirichlet_condition_brick :
public virtual_brick {
4303 bool normal_component;
4304 const mesh_fem *mf_mult_;
4310 virtual void asm_real_tangent_terms(
const model &md,
size_type ib,
4311 const model::varnamelist &vl,
4312 const model::varnamelist &dl,
4313 const model::mimlist &mims,
4314 model::real_matlist &matl,
4315 model::real_veclist &vecl,
4316 model::real_veclist &,
4318 build_version version)
const {
4319 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
4320 "Dirichlet condition brick has one and only one term");
4321 GMM_ASSERT1(mims.size() == 1,
4322 "Dirichlet condition brick need one and only one mesh_im");
4323 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() <= 3,
4324 "Wrong number of variables for Dirichlet condition brick");
4326 model_real_sparse_matrix& rB = rB_th;
4327 model_real_plain_vector& rV = rV_th;
4329 bool penalized = (vl.size() == 1);
4330 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4331 const mesh_fem &mf_mult = penalized ? (mf_mult_ ? *mf_mult_ : mf_u)
4332 : md.mesh_fem_of_variable(vl[1]);
4333 const mesh_im &mim = *mims[0];
4334 const model_real_plain_vector *A = 0, *COEFF = 0, *H = 0;
4335 const mesh_fem *mf_data = 0, *mf_H = 0;
4336 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
4337 || (penalized && md.is_var_newer_than_brick(dl[0], ib));
4340 COEFF = &(md.real_variable(dl[0]));
4341 GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
4342 "Data for coefficient should be a scalar");
4345 size_type s = 0, ind = (penalized ? 1 : 0);
4346 if (dl.size() > ind) {
4347 A = &(md.real_variable(dl[ind]));
4348 mf_data = md.pmesh_fem_of_variable(dl[ind]);
4349 s = gmm::vect_size(*A);
4350 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4351 size_type ss = ((normal_component) ? 1 : mf_u.get_qdim());
4352 GMM_ASSERT1(s == ss, dl[ind] <<
": bad format of "
4353 "Dirichlet data. Detected dimension is " << s
4354 <<
" should be " << ss);
4357 if (dl.size() > ind + 1) {
4358 GMM_ASSERT1(H_version,
4359 "Wrong number of data for Dirichlet condition brick");
4360 H = &(md.real_variable(dl[ind+1]));
4361 mf_H = md.pmesh_fem_of_variable(dl[ind+1]);
4362 s = gmm::vect_size(*A);
4364 s = s * mf_H->get_qdim() / mf_H->nb_dof();
4365 GMM_ASSERT1(mf_H->get_qdim() == 1,
"Implemented only for mf_H "
4366 "a scalar finite element method");
4368 GMM_ASSERT1(s = gmm::sqr(mf_u.get_qdim()),
4369 dl[ind+1] <<
": bad format of Dirichlet data. "
4370 "Detected dimension is " << s <<
" should be "
4371 <<
size_type(gmm::sqr(mf_u.get_qdim())));
4374 mesh_region rg(region);
4375 mim.linked_mesh().intersect_with_mpi_region(rg);
4377 if (recompute_matrix) {
4378 model_real_sparse_matrix *B = &(matl[0]);
4379 if (penalized && (&mf_mult != &mf_u)) {
4386 GMM_TRACE2(
"Mass term assembly for Dirichlet condition");
4387 if (H_version || normal_component) {
4388 ga_workspace workspace;
4389 gmm::sub_interval Imult(0, mf_mult.nb_dof()), Iu(0, mf_u.nb_dof());
4390 base_vector u(mf_u.nb_dof());
4391 base_vector mult(mf_mult.nb_dof());
4392 workspace.add_fem_variable(
"u", mf_u, Iu, u);
4393 workspace.add_fem_variable(
"mult", mf_mult, Imult, mult);
4394 auto expression = std::string{};
4396 if (mf_H) workspace.add_fem_constant(
"A", *mf_H, *H);
4397 else workspace.add_fixed_size_constant(
"A", *H);
4398 expression = (mf_u.get_qdim() == 1) ?
4399 "Test_mult . (A . Test2_u)"
4401 "Test_mult. (Reshape(A, qdim(u), qdim(u)) . Test2_u)";
4402 }
else if (normal_component) {
4403 expression =
"Test_mult . (Test2_u . Normal)";
4405 workspace.add_expression(expression, mim, rg);
4406 workspace.set_assembled_matrix(*B);
4407 workspace.assembly(2);
4412 if (penalized && (&mf_mult != &mf_u)) {
4413 GMM_ASSERT1(MPI_IS_MASTER(),
"Sorry, the penalized option "
4414 "filtered by a multiplier space is not parallelized");
4415 gmm::mult(gmm::transposed(rB), rB, matl[0]);
4416 gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
4417 }
else if (penalized) {
4418 gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
4422 if (dl.size() > ind) {
4423 GMM_TRACE2(
"Source term assembly for Dirichlet condition");
4425 if (penalized && (&mf_mult != &mf_u)) {
4439 if (penalized && (&mf_mult != &mf_u)) {
4440 gmm::mult(gmm::transposed(rB), rV, vecl[0]);
4441 gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
4442 rV = model_real_plain_vector();
4443 }
else if (penalized)
4444 gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
4449 void real_post_assembly_in_serial(
const model &md,
size_type ib,
4450 const model::varnamelist &,
4451 const model::varnamelist &,
4452 const model::mimlist &,
4453 model::real_matlist &,
4454 model::real_veclist &vecl,
4455 model::real_veclist &,
4457 build_version)
const override {
4458 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4461 virtual void asm_complex_tangent_terms(
const model &md,
size_type ib,
4462 const model::varnamelist &vl,
4463 const model::varnamelist &dl,
4464 const model::mimlist &mims,
4465 model::complex_matlist &matl,
4466 model::complex_veclist &vecl,
4467 model::complex_veclist &,
4469 build_version version)
const {
4470 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
4471 "Dirichlet condition brick has one and only one term");
4472 GMM_ASSERT1(mims.size() == 1,
4473 "Dirichlet condition brick need one and only one mesh_im");
4474 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() <= 3,
4475 "Wrong number of variables for Dirichlet condition brick");
4477 model_complex_sparse_matrix& cB = cB_th;
4478 model_complex_plain_vector& cV = cV_th;
4480 bool penalized = (vl.size() == 1);
4481 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4482 const mesh_fem &mf_mult = penalized ? (mf_mult_ ? *mf_mult_ : mf_u)
4483 : md.mesh_fem_of_variable(vl[1]);
4484 const mesh_im &mim = *mims[0];
4485 const model_complex_plain_vector *A = 0, *COEFF = 0, *H = 0;
4486 const mesh_fem *mf_data = 0, *mf_H = 0;
4487 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
4488 || (penalized && md.is_var_newer_than_brick(dl[0], ib));
4491 COEFF = &(md.complex_variable(dl[0]));
4492 GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
4493 "Data for coefficient should be a scalar");
4496 size_type s = 0, ind = (penalized ? 1 : 0);
4497 if (dl.size() > ind) {
4498 A = &(md.complex_variable(dl[ind]));
4499 mf_data = md.pmesh_fem_of_variable(dl[ind]);
4500 s = gmm::vect_size(*A);
4501 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4502 size_type ss = s * ((normal_component) ? mf_u.linked_mesh().dim() : 1);
4503 GMM_ASSERT1(mf_u.get_qdim() == ss,
4504 dl[ind] <<
": bad format of Dirichlet data. "
4505 "Detected dimension is " << ss <<
" should be "
4509 if (dl.size() > ind + 1) {
4510 GMM_ASSERT1(H_version,
4511 "Wrong number of data for Dirichlet condition brick");
4512 H = &(md.complex_variable(dl[ind+1]));
4513 mf_H = md.pmesh_fem_of_variable(dl[ind+1]);
4514 s = gmm::vect_size(*A);
4516 s = s * mf_H->get_qdim() / mf_H->nb_dof();
4517 GMM_ASSERT1(mf_H->get_qdim() == 1,
"Implemented only for mf_H "
4518 "a scalar finite element method");
4520 GMM_ASSERT1(s = gmm::sqr(mf_u.get_qdim()),
4521 dl[ind+1] <<
": bad format of Dirichlet data. "
4522 "Detected dimension is " << s <<
" should be "
4523 <<
size_type(gmm::sqr(mf_u.get_qdim())));
4526 mesh_region rg(region);
4527 mim.linked_mesh().intersect_with_mpi_region(rg);
4529 if (recompute_matrix) {
4530 model_complex_sparse_matrix *B = &(matl[0]);
4531 if (penalized && (&mf_mult != &mf_u)) {
4538 GMM_TRACE2(
"Mass term assembly for Dirichlet condition");
4540 if (mf_u.get_qdim() == 1)
4541 asm_real_or_complex_1_param_mat(*B, mim, mf_mult, mf_H, *H, rg,
4542 "(A*Test_u).Test2_u");
4544 asm_real_or_complex_1_param_mat(*B, mim, mf_mult, mf_H, *H, rg,
4545 "(Reshape(A,qdim(u),qdim(u))*Test2_u).Test_u");
4561 else if (normal_component) {
4562 ga_workspace workspace;
4563 gmm::sub_interval Imult(0, mf_mult.nb_dof()), Iu(0, mf_u.nb_dof());
4564 base_vector mult(mf_mult.nb_dof()), u(mf_u.nb_dof());
4565 workspace.add_fem_variable(
"mult", mf_mult, Imult, mult);
4566 workspace.add_fem_variable(
"u", mf_u, Iu, u);
4567 workspace.add_expression(
"Test_mult.(Test2_u.Normal)", mim, rg);
4568 model_real_sparse_matrix BB(mf_mult.nb_dof(), mf_u.nb_dof());
4569 workspace.set_assembled_matrix(BB);
4570 workspace.assembly(2);
4586 if (penalized && (&mf_mult != &mf_u)) {
4587 gmm::mult(gmm::transposed(cB), cB, matl[0]);
4588 gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
4589 }
else if (penalized) {
4590 gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
4594 if (dl.size() > ind) {
4595 GMM_TRACE2(
"Source term assembly for Dirichlet condition");
4597 if (penalized && (&mf_mult != &mf_u)) {
4611 if (penalized && (&mf_mult != &mf_u)) {
4612 gmm::mult(gmm::transposed(cB), cV, vecl[0]);
4613 gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
4614 cV = model_complex_plain_vector();
4615 }
else if (penalized)
4616 gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
4620 void complex_post_assembly_in_serial(
const model &md,
size_type ib,
4621 const model::varnamelist &,
4622 const model::varnamelist &,
4623 const model::mimlist &,
4624 model::complex_matlist &,
4625 model::complex_veclist &vecl,
4626 model::complex_veclist &,
4628 build_version)
const override {
4629 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4633 virtual std::string declare_volume_assembly_string
4634 (
const model &,
size_type,
const model::varnamelist &,
4635 const model::varnamelist &)
const {
4636 return std::string();
4639 Dirichlet_condition_brick(
bool penalized,
bool H_version_,
4640 bool normal_component_,
4641 const mesh_fem *mf_mult__ = 0) {
4642 mf_mult_ = mf_mult__;
4643 H_version = H_version_;
4644 normal_component = normal_component_;
4645 GMM_ASSERT1(!(H_version && normal_component),
"Bad Dirichlet version");
4646 set_flags(penalized ?
"Dirichlet with penalization brick"
4647 :
"Dirichlet with multipliers brick",
4657 const std::string &multname,
size_type region,
4658 const std::string &dataname) {
4659 pbrick pbr = std::make_shared<Dirichlet_condition_brick>(
false,
false,
false);
4661 tl.push_back(model::term_description(multname, varname,
true));
4662 model::varnamelist vl(1, varname);
4663 vl.push_back(multname);
4664 model::varnamelist dl;
4665 if (dataname.size()) dl.push_back(dataname);
4666 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4672 const std::string &dataname) {
4673 std::string multname = md.
new_name(
"mult_on_" + varname);
4676 (md, mim, varname, multname, region, dataname);
4682 const std::string &dataname) {
4687 (md, mim, varname, mf_mult, region, dataname);
4696 scalar_type penalisation_coeff,
size_type region,
4697 const std::string &dataname,
const mesh_fem *mf_mult) {
4698 std::string coeffname = md.
new_name(
"penalization_on_" + varname);
4704 pbrick pbr = std::make_shared<Dirichlet_condition_brick>
4705 (
true,
false,
false, mf_mult);
4707 tl.push_back(model::term_description(varname, varname,
true));
4708 model::varnamelist vl(1, varname);
4709 model::varnamelist dl(1, coeffname);
4710 if (dataname.size()) dl.push_back(dataname);
4711 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4716 const std::string &multname,
size_type region,
4717 const std::string &dataname) {
4718 pbrick pbr = std::make_shared<Dirichlet_condition_brick>(
false,
false,
true);
4720 tl.push_back(model::term_description(multname, varname,
true));
4721 model::varnamelist vl(1, varname);
4722 vl.push_back(multname);
4723 model::varnamelist dl;
4724 if (dataname.size()) dl.push_back(dataname);
4725 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4731 const std::string &dataname) {
4732 std::string multname = md.
new_name(
"mult_on_" + varname);
4735 (md, mim, varname, multname, region, dataname);
4741 const std::string &dataname) {
4745 (md, mim, varname, mf_mult, region, dataname);
4750 scalar_type penalisation_coeff,
size_type region,
4751 const std::string &dataname,
const mesh_fem *mf_mult) {
4752 std::string coeffname = md.
new_name(
"penalization_on_" + varname);
4758 pbrick pbr = std::make_shared<Dirichlet_condition_brick>
4759 (
true,
false,
true, mf_mult);
4761 tl.push_back(model::term_description(varname, varname,
true));
4762 model::varnamelist vl(1, varname);
4763 model::varnamelist dl(1, coeffname);
4764 if (dataname.size()) dl.push_back(dataname);
4765 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4771 const std::string &multname,
size_type region,
4772 const std::string &dataname,
const std::string &Hname) {
4773 pbrick pbr = std::make_shared<Dirichlet_condition_brick>(
false,
true,
false);
4775 tl.push_back(model::term_description(multname, varname,
true));
4776 model::varnamelist vl(1, varname);
4777 vl.push_back(multname);
4778 model::varnamelist dl;
4779 dl.push_back(dataname);
4780 dl.push_back(Hname);
4781 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4787 const std::string &dataname,
const std::string &Hname) {
4788 std::string multname = md.
new_name(
"mult_on_" + varname);
4791 (md, mim, varname, multname, region, dataname, Hname);
4797 const std::string &dataname,
const std::string &Hname) {
4802 (md, mim, varname, mf_mult, region, dataname, Hname);
4807 scalar_type penalisation_coeff,
size_type region,
4808 const std::string &dataname,
const std::string &Hname,
4810 std::string coeffname = md.
new_name(
"penalization_on_" + varname);
4816 pbrick pbr = std::make_shared<Dirichlet_condition_brick>
4817 (
true,
true,
false, mf_mult);
4819 tl.push_back(model::term_description(varname, varname,
true));
4820 model::varnamelist vl(1, varname);
4821 model::varnamelist dl(1, coeffname);
4822 dl.push_back(dataname);
4823 dl.push_back(Hname);
4824 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4828 scalar_type penalisation_coeff) {
4832 GMM_ASSERT1(gmm::vect_size(d)==1,
4833 "Wrong coefficient size, may be not a Dirichlet brick "
4834 "with penalization");
4835 d[0] = penalisation_coeff;
4839 GMM_ASSERT1(gmm::vect_size(d)==1,
4840 "Wrong coefficient size, may be not a Dirichlet brick "
4841 "with penalization");
4842 d[0] = penalisation_coeff;
4852 struct simplification_Dirichlet_condition_brick :
public virtual_brick {
4854 virtual void asm_real_tangent_terms(
const model &md,
size_type ,
4855 const model::varnamelist &vl,
4856 const model::varnamelist &dl,
4857 const model::mimlist &mims,
4858 model::real_matlist &matl,
4859 model::real_veclist &vecl,
4860 model::real_veclist &,
4862 build_version )
const {
4863 if (MPI_IS_MASTER()) {
4865 GMM_ASSERT1(vecl.size() == 0 && matl.size() == 0,
4866 "Dirichlet condition brick by simplification has no term");
4867 GMM_ASSERT1(mims.size() == 0,
4868 "Dirichlet condition brick by simplification need no "
4870 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
4871 "Wrong number of variables for Dirichlet condition brick "
4872 "by simplification");
4874 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4875 const model_real_plain_vector *A = 0;
4876 const mesh_fem *mf_data = 0;
4879 if (dl.size() == 1) {
4880 A = &(md.real_variable(dl[0]));
4881 mf_data = md.pmesh_fem_of_variable(dl[0]);
4884 GMM_ASSERT1(mf_data == &mf_u,
"Sorry, for this brick, the data has"
4885 " to be defined on the same f.e.m. as the unknown");
4887 s = gmm::vect_size(*A);
4888 GMM_ASSERT1(mf_u.get_qdim() == s,
": bad format of "
4889 "Dirichlet data. Detected dimension is " << s
4890 <<
" should be " <<
size_type(mf_u.get_qdim()));
4894 mesh_region rg(region);
4897 if (mf_u.get_qdim() > 1 || (!mf_data && A)) {
4898 for (mr_visitor i(rg, mf_u.linked_mesh()); !i.finished(); ++i) {
4899 pfem pf = mf_u.fem_of_element(i.cv());
4901 GMM_ASSERT1(pf->target_dim() == 1,
4902 "Intrinsically vectorial fems are not allowed");
4903 GMM_ASSERT1(mf_data || pf->is_lagrange(),
"Constant Dirichlet "
4904 "data allowed for lagrange fems only");
4909 dal::bit_vector dofs = mf_u.dof_on_region(rg);
4911 if (A && !mf_data) {
4912 GMM_ASSERT1(dofs.card() % s == 0,
"Problem with dof vectorization");
4915 for (dal::bv_visitor i(dofs); !i.finished(); ++i) {
4917 if (A) val = (mf_data ? (*A)[i] : (*A)[i%s]);
4918 md.add_real_dof_constraint(vl[0], i, val);
4923 virtual void asm_complex_tangent_terms(
const model &md,
size_type ,
4924 const model::varnamelist &vl,
4925 const model::varnamelist &dl,
4926 const model::mimlist &mims,
4927 model::complex_matlist &matl,
4928 model::complex_veclist &vecl,
4929 model::complex_veclist &,
4931 build_version )
const {
4932 if (MPI_IS_MASTER()) {
4933 GMM_ASSERT1(vecl.size() == 0 && matl.size() == 0,
4934 "Dirichlet condition brick by simplification has no term");
4935 GMM_ASSERT1(mims.size() == 0,
4936 "Dirichlet condition brick by simplification need no "
4938 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
4939 "Wrong number of variables for Dirichlet condition brick "
4940 "by simplification");
4942 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4943 const model_complex_plain_vector *A = 0;
4944 const mesh_fem *mf_data = 0;
4947 if (dl.size() == 1) {
4948 A = &(md.complex_variable(dl[0]));
4949 mf_data = md.pmesh_fem_of_variable(dl[0]);
4952 GMM_ASSERT1(mf_data == &mf_u,
"Sorry, for this brick, the data has"
4953 " to be define on the same f.e.m. than the unknown");
4955 s = gmm::vect_size(*A);
4956 GMM_ASSERT1(mf_u.get_qdim() == s,
": bad format of "
4957 "Dirichlet data. Detected dimension is " << s
4958 <<
" should be " <<
size_type(mf_u.get_qdim()));
4962 mesh_region rg(region);
4965 if (mf_u.get_qdim() > 1 || (!mf_data && A)) {
4966 for (mr_visitor i(rg, mf_u.linked_mesh()); !i.finished(); ++i) {
4967 pfem pf = mf_u.fem_of_element(i.cv());
4969 GMM_ASSERT1(pf->target_dim() == 1,
4970 "Intrinsically vectorial fems are not allowed");
4971 GMM_ASSERT1(mf_data || pf->is_lagrange(),
"Constant Dirichlet "
4972 "data allowed for lagrange fems only");
4977 dal::bit_vector dofs = mf_u.dof_on_region(rg);
4979 if (A && !mf_data) {
4980 GMM_ASSERT1(dofs.card() % s == 0,
"Problem with dof vectorization");
4983 for (dal::bv_visitor i(dofs); !i.finished(); ++i) {
4984 complex_type val(0);
4985 if (A) val = (mf_data ? (*A)[i] : (*A)[i%s]);
4986 md.add_complex_dof_constraint(vl[0], i, val);
4992 virtual std::string declare_volume_assembly_string
4993 (
const model &,
size_type,
const model::varnamelist &,
4994 const model::varnamelist &)
const {
4995 return std::string();
4998 simplification_Dirichlet_condition_brick() {
4999 set_flags(
"Dirichlet with simplification brick",
5009 size_type region,
const std::string &dataname) {
5010 pbrick pbr = std::make_shared<simplification_Dirichlet_condition_brick>();
5012 model::varnamelist vl(1, varname);
5013 model::varnamelist dl;
5014 if (dataname.size()) dl.push_back(dataname);
5015 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(), region);
5026 const std::string &Neumannterm,
5027 const std::string &datagamma0,
size_type region, scalar_type theta_,
5028 const std::string &datag) {
5029 std::string theta = std::to_string(theta_);
5031 ga_workspace workspace(md, ga_workspace::inherit::ALL);
5032 size_type order = workspace.add_expression(Neumannterm, mim, region, 1);
5033 GMM_ASSERT1(order == 0,
"Wrong expression of the Neumann term");
5034 bool is_lin = workspace.is_linear(1);
5036 std::string condition =
"("+varname + (datag.size() ?
"-("+datag+
"))":
")");
5037 std::string gamma =
"(("+datagamma0+
")*element_size)";
5038 std::string r =
"(1/"+gamma+
")";
5039 std::string expr =
"("+r+
"*"+condition+
"-("+Neumannterm+
")).Test_"+varname;
5040 if (theta_ != scalar_type(0)) {
5041 std::string derivative_Neumann = workspace.extract_order1_term(varname);
5042 if (derivative_Neumann.size())
5043 expr+=
"-"+theta+
"*"+condition+
".("+derivative_Neumann+
")";
5051 "Dirichlet condition with Nitsche's method");
5054 "Dirichlet condition with Nitsche's method");
5060 const std::string &Neumannterm,
5061 const std::string &datagamma0,
size_type region, scalar_type theta_,
5062 const std::string &datag) {
5063 std::string theta = std::to_string(theta_);
5065 ga_workspace workspace(md, ga_workspace::inherit::ALL);
5066 size_type order = workspace.add_expression(Neumannterm, mim, region, 1);
5067 GMM_ASSERT1(order == 0,
"Wrong expression of the Neumann term");
5068 bool is_lin = workspace.is_linear(1);
5070 std::string condition =
"("+varname+
".Normal"
5071 + (datag.size() ?
"-("+datag+
"))":
")");
5072 std::string gamma =
"(("+datagamma0+
")*element_size)";
5073 std::string r =
"(1/"+gamma+
")";
5074 std::string expr =
"("+r+
"*"+condition+
"-Normal.("+Neumannterm
5075 +
"))*(Normal.Test_"+varname+
")";
5076 if (theta_ != scalar_type(0)) {
5077 std::string derivative_Neumann = workspace.extract_order1_term(varname);
5078 if (derivative_Neumann.size())
5079 expr+=
"-"+theta+
"*"+condition+
"*Normal.("
5080 +derivative_Neumann+
")";
5084 "Dirichlet condition with Nitsche's method");
5087 "Dirichlet condition with Nitsche's method");
5093 const std::string &Neumannterm,
5094 const std::string &datagamma0,
size_type region, scalar_type theta_,
5095 const std::string &datag,
const std::string &dataH) {
5096 std::string theta = std::to_string(theta_);
5098 ga_workspace workspace(md, ga_workspace::inherit::ALL);
5099 size_type order = workspace.add_expression(Neumannterm, mim, region, 1);
5100 GMM_ASSERT1(order == 0,
"Wrong expression of the Neumann term");
5101 bool is_lin = workspace.is_linear(1);
5103 std::string condition =
"(("+dataH+
")*"+varname
5104 + (datag.size() ?
"-("+datag+
"))":
")");
5105 std::string gamma =
"(("+datagamma0+
")*element_size)";
5106 std::string r =
"(1/"+gamma+
")";
5107 std::string expr =
"("+r+
"*"+condition+
"-("+dataH+
")*("+Neumannterm
5108 +
"))*(("+dataH+
")*Test_"+varname+
")";
5109 if (theta_ != scalar_type(0)) {
5110 std::string derivative_Neumann = workspace.extract_order1_term(varname);
5111 if (derivative_Neumann.size())
5112 expr+=
"-"+theta+
"*"+condition+
"*(("+dataH+
")*("
5113 +derivative_Neumann+
"))";
5117 "Dirichlet condition with Nitsche's method");
5120 "Dirichlet condition with Nitsche's method");
5132 struct pointwise_constraints_brick :
public virtual_brick {
5134 mutable gmm::row_matrix<model_real_sparse_vector> rB;
5135 mutable gmm::row_matrix<model_complex_sparse_vector> cB;
5137 virtual void real_pre_assembly_in_serial(
const model &md,
size_type ib,
5138 const model::varnamelist &vl,
5139 const model::varnamelist &dl,
5140 const model::mimlist &mims,
5141 model::real_matlist &matl,
5142 model::real_veclist &vecl,
5143 model::real_veclist &,
5145 build_version version)
const {
5146 if (MPI_IS_MASTER()) {
5148 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5149 "Pointwize constraints brick has only one term");
5150 GMM_ASSERT1(mims.size() == 0,
5151 "Pointwize constraints brick does not need a mesh_im");
5152 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2,
5153 "Wrong number of variables for pointwize constraints brick");
5154 bool penalized = (vl.size() == 1);
5155 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5156 dim_type N = mf_u.linked_mesh().dim(), Q = mf_u.get_qdim(), ind_pt = 0;
5158 GMM_ASSERT1(dl.size() == dlsize || dl.size() == dlsize+1,
5159 "Wrong number of data for pointwize constraints brick");
5162 const model_real_plain_vector *COEFF = 0;
5164 COEFF = &(md.real_variable(dl[0]));
5166 GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
5167 "Data for coefficient should be a scalar");
5170 const model_real_plain_vector &PT = md.real_variable(dl[ind_pt]);
5171 size_type nb_co = gmm::vect_size(PT) / N;
5173 dim_type ind_unitv = dim_type((Q > 1) ? ind_pt+1 : 0);
5174 const model_real_plain_vector &unitv =md.real_variable(dl[ind_unitv]);
5175 GMM_ASSERT1((!ind_unitv || gmm::vect_size(unitv) == nb_co * Q),
5176 "Wrong size for vector of unit vectors");
5178 dim_type ind_rhs = dim_type((Q > 1) ? ind_pt+2 : ind_pt+1);
5179 if (dl.size() <
size_type(ind_rhs + 1)) ind_rhs = 0;
5180 const model_real_plain_vector &rhs = md.real_variable(dl[ind_rhs]);
5181 GMM_ASSERT1((!ind_rhs || gmm::vect_size(rhs) == nb_co),
5182 "Wrong size for vector of rhs");
5184 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
5185 || (penalized && (md.is_var_newer_than_brick(dl[ind_pt], ib)
5186 || md.is_var_newer_than_brick(dl[ind_unitv], ib)
5187 || md.is_var_newer_than_brick(dl[ind_rhs], ib)));
5189 if (recompute_matrix) {
5190 gmm::row_matrix<model_real_sparse_vector> BB(nb_co*Q, mf_u.nb_dof());
5193 dal::bit_vector dof_untouched;
5194 getfem::mesh_trans_inv mti(mf_u.linked_mesh());
5197 gmm::copy(gmm::sub_vector(PT, gmm::sub_interval(i*N, N)), pt);
5200 gmm::row_matrix<model_real_sparse_vector> &BBB = ((Q > 1) ? BB : rB);
5201 model_real_plain_vector vv;
5202 interpolation(mf_u, mti, vv, vv, BBB, 1, 1, &dof_untouched);
5203 GMM_ASSERT1(dof_untouched.card() == 0,
5204 "Pointwize constraints : some of the points are outside "
5205 "the mesh: " << dof_untouched);
5210 gmm::add(gmm::scaled(gmm::mat_row(BB, i*Q+q), unitv[i*Q+q]),
5211 gmm::mat_row(rB, i));
5214 gmm::mult(gmm::transposed(rB), rB, matl[0]);
5215 gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
5217 gmm::copy(rB, matl[0]);
5222 gmm::mult(gmm::transposed(rB), rhs, vecl[0]);
5223 gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
5225 else gmm::copy(rhs, vecl[0]);
5231 virtual void complex_pre_assembly_in_serial(
const model &md,
size_type ib,
5232 const model::varnamelist &vl,
5233 const model::varnamelist &dl,
5234 const model::mimlist &mims,
5235 model::complex_matlist &matl,
5236 model::complex_veclist &vecl,
5237 model::complex_veclist &,
5239 build_version version)
const {
5240 if (MPI_IS_MASTER()) {
5241 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5242 "Pointwize constraints brick only one term");
5243 GMM_ASSERT1(mims.size() == 0,
5244 "Pointwize constraints brick does not need a mesh_im");
5245 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2,
5246 "Wrong number of variables for pointwize constraints brick");
5247 bool penalized = (vl.size() == 1);
5248 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5249 dim_type N = mf_u.linked_mesh().dim(), Q = mf_u.get_qdim(), ind_pt = 0;
5251 GMM_ASSERT1(dl.size() == dlsize || dl.size() == dlsize+1,
5252 "Wrong number of data for pointwize constraints brick");
5255 const model_complex_plain_vector *COEFF = 0;
5257 COEFF = &(md.complex_variable(dl[0]));
5259 GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
5260 "Data for coefficient should be a scalar");
5263 const model_complex_plain_vector &PT = md.complex_variable(dl[ind_pt]);
5264 size_type nb_co = gmm::vect_size(PT) / N;
5266 dim_type ind_unitv = dim_type((Q > 1) ? ind_pt+1 : 0);
5267 const model_complex_plain_vector &unitv
5268 = md.complex_variable(dl[ind_unitv]);
5269 GMM_ASSERT1((!ind_unitv || gmm::vect_size(unitv) == nb_co * Q),
5270 "Wrong size for vector of unit vectors");
5272 dim_type ind_rhs = dim_type((Q > 1) ? ind_pt+2 : ind_pt+1);
5273 if (dl.size() <
size_type(ind_rhs + 1)) ind_rhs = 0;
5274 const model_complex_plain_vector &rhs
5275 = md.complex_variable(dl[ind_rhs]);
5276 GMM_ASSERT1((!ind_rhs || gmm::vect_size(rhs) == nb_co),
5277 "Wrong size for vector of rhs");
5279 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
5280 || (penalized && (md.is_var_newer_than_brick(dl[ind_pt], ib)
5281 || md.is_var_newer_than_brick(dl[ind_unitv], ib)
5282 || md.is_var_newer_than_brick(dl[ind_rhs], ib)));
5284 if (recompute_matrix) {
5285 gmm::row_matrix<model_complex_sparse_vector>
5286 BB(nb_co*Q,mf_u.nb_dof());
5288 dal::bit_vector dof_untouched;
5289 getfem::mesh_trans_inv mti(mf_u.linked_mesh());
5292 gmm::copy(gmm::real_part
5293 (gmm::sub_vector(PT, gmm::sub_interval(i*N, N))), pt);
5296 gmm::row_matrix<model_complex_sparse_vector> &BBB = ((Q > 1) ? BB :cB);
5297 model_complex_plain_vector vv;
5298 interpolation(mf_u, mti, vv, vv, BBB, 1, 1, &dof_untouched);
5299 GMM_ASSERT1(dof_untouched.card() == 0,
5300 "Pointwize constraints : some of the points are outside "
5301 "the mesh: " << dof_untouched);
5306 gmm::add(gmm::scaled(gmm::mat_row(BB, i*Q+q), unitv[i*Q+q]),
5307 gmm::mat_row(cB, i));
5311 gmm::mult(gmm::transposed(cB), cB, matl[0]);
5312 gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
5314 gmm::copy(cB, matl[0]);
5320 gmm::mult(gmm::transposed(cB), rhs, vecl[0]);
5321 gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
5323 else gmm::copy(rhs, vecl[0]);
5329 virtual std::string declare_volume_assembly_string
5330 (
const model &,
size_type,
const model::varnamelist &,
5331 const model::varnamelist &)
const {
5332 return std::string();
5335 pointwise_constraints_brick(
bool penalized) {
5336 set_flags(penalized ?
"Pointwise cosntraints with penalization brick"
5337 :
"Pointwise cosntraints with multipliers brick",
5348 scalar_type penalisation_coeff,
const std::string &dataname_pt,
5349 const std::string &dataname_unitv,
const std::string &dataname_val) {
5350 std::string coeffname = md.
new_name(
"penalization_on_" + varname);
5356 pbrick pbr = std::make_shared<pointwise_constraints_brick>(
true);
5358 tl.push_back(model::term_description(varname, varname,
true));
5359 model::varnamelist vl(1, varname);
5360 model::varnamelist dl(1, coeffname);
5361 dl.push_back(dataname_pt);
5363 if (mf_u.
get_qdim() > 1) dl.push_back(dataname_unitv);
5364 if (dataname_val.size() > 0) dl.push_back(dataname_val);
5370 const std::string &multname,
const std::string &dataname_pt,
5371 const std::string &dataname_unitv,
const std::string &dataname_val) {
5372 pbrick pbr = std::make_shared<pointwise_constraints_brick>(
false);
5374 tl.push_back(model::term_description(multname, varname,
true));
5375 model::varnamelist vl(1, varname);
5376 vl.push_back(multname);
5377 model::varnamelist dl(1, dataname_pt);
5379 if (mf_u.
get_qdim() > 1) dl.push_back(dataname_unitv);
5380 if (dataname_val.size() > 0) dl.push_back(dataname_val);
5385 (
model &md,
const std::string &varname,
const std::string &dataname_pt,
5386 const std::string &dataname_unitv,
const std::string &dataname_val) {
5387 std::string multname = md.
new_name(
"mult_on_" + varname);
5395 (md, varname, multname, dataname_pt, dataname_unitv, dataname_val);
5405 struct Helmholtz_brick :
public virtual_brick {
5407 virtual void asm_real_tangent_terms(
const model &md,
size_type,
5408 const model::varnamelist &vl,
5409 const model::varnamelist &dl,
5410 const model::mimlist &mims,
5411 model::real_matlist &matl,
5412 model::real_veclist &,
5413 model::real_veclist &,
5415 build_version)
const {
5416 GMM_ASSERT1(matl.size() == 1,
5417 "Helmholtz brick has one and only one term");
5418 GMM_ASSERT1(mims.size() == 1,
5419 "Helmholtz brick need one and only one mesh_im");
5420 GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
5421 "Wrong number of variables for Helmholtz brick");
5423 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5424 const mesh &m = mf_u.linked_mesh();
5426 GMM_ASSERT1(Q == 1,
"Helmholtz brick is only for scalar field, sorry.");
5427 const mesh_im &mim = *mims[0];
5428 const mesh_fem *mf_a = 0;
5429 mesh_region rg(region);
5430 m.intersect_with_mpi_region(rg);
5431 const model_real_plain_vector *A = &(md.real_variable(dl[0]));
5432 mf_a = md.pmesh_fem_of_variable(dl[0]);
5433 s = gmm::vect_size(*A);
5434 if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
5437 GMM_TRACE2(
"Stiffness matrix assembly for Helmholtz problem");
5438 gmm::clear(matl[0]);
5439 model_real_plain_vector A2(gmm::vect_size(*A));
5440 for (
size_type i=0; i < gmm::vect_size(*A); ++i)
5441 A2[i] = gmm::sqr((*A)[i]);
5447 GMM_ASSERT1(
false,
"Bad format Helmholtz brick coefficient");
5450 virtual void asm_complex_tangent_terms(
const model &md,
size_type,
5451 const model::varnamelist &vl,
5452 const model::varnamelist &dl,
5453 const model::mimlist &mims,
5454 model::complex_matlist &matl,
5455 model::complex_veclist &,
5456 model::complex_veclist &,
5458 build_version)
const {
5459 GMM_ASSERT1(matl.size() == 1,
5460 "Helmholtz brick has one and only one term");
5461 GMM_ASSERT1(mims.size() == 1,
5462 "Helmholtz brick need one and only one mesh_im");
5463 GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
5464 "Wrong number of variables for Helmholtz brick");
5466 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5467 const mesh &m = mf_u.linked_mesh();
5469 GMM_ASSERT1(Q == 1,
"Helmholtz brick is only for scalar field, sorry.");
5470 const mesh_im &mim = *mims[0];
5471 const mesh_fem *mf_a = 0;
5472 mesh_region rg(region);
5473 m.intersect_with_mpi_region(rg);
5474 const model_complex_plain_vector *A = &(md.complex_variable(dl[0]));
5475 mf_a = md.pmesh_fem_of_variable(dl[0]);
5476 s = gmm::vect_size(*A);
5477 if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
5480 GMM_TRACE2(
"Stiffness matrix assembly for Helmholtz problem");
5482 model_complex_plain_vector A2(gmm::vect_size(*A));
5483 for (
size_type i=0; i < gmm::vect_size(*A); ++i)
5484 A2[i] = gmm::sqr((*A)[i]);
5490 GMM_ASSERT1(
false,
"Bad format Helmholtz brick coefficient");
5494 set_flags(
"Helmholtz",
true ,
5502 const std::string &varname,
5503 const std::string &dataexpr,
5506 pbrick pbr = std::make_shared<Helmholtz_brick>();
5508 tl.push_back(model::term_description(varname, varname,
true));
5509 return md.
add_brick(pbr, model::varnamelist(1, varname),
5510 model::varnamelist(1, dataexpr), tl,
5511 model::mimlist(1, &mim), region);
5513 std::string test_varname
5514 =
"Test_" + sup_previous_and_dot_to_varname(varname);
5515 std::string expr =
"Grad_"+varname+
".Grad_"+test_varname
5516 +
" + sqr("+dataexpr+
")*"+varname+
"*"+test_varname;
5522 "Helmholtz (nonlinear)");
5535 struct Fourier_Robin_brick :
public virtual_brick {
5537 virtual void asm_real_tangent_terms(
const model &md,
size_type,
5538 const model::varnamelist &vl,
5539 const model::varnamelist &dl,
5540 const model::mimlist &mims,
5541 model::real_matlist &matl,
5542 model::real_veclist &,
5543 model::real_veclist &,
5545 build_version)
const {
5546 GMM_ASSERT1(matl.size() == 1,
5547 "Fourier-Robin brick has one and only one term");
5548 GMM_ASSERT1(mims.size() == 1,
5549 "Fourier-Robin brick need one and only one mesh_im");
5550 GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
5551 "Wrong number of variables for Fourier-Robin brick");
5553 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5554 const mesh &m = mf_u.linked_mesh();
5556 const mesh_im &mim = *mims[0];
5557 const mesh_fem *mf_a = 0;
5558 mesh_region rg(region);
5559 m.intersect_with_mpi_region(rg);
5560 const model_real_plain_vector *A = &(md.real_variable(dl[0]));
5561 mf_a = md.pmesh_fem_of_variable(dl[0]);
5562 s = gmm::vect_size(*A);
5563 if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
5564 GMM_ASSERT1(s == Q*Q,
5565 "Bad format Fourier-Robin brick coefficient");
5567 GMM_TRACE2(
"Fourier-Robin term assembly");
5568 gmm::clear(matl[0]);
5572 asm_homogeneous_qu_term(matl[0], mim, mf_u, *A, rg);
5575 virtual void asm_complex_tangent_terms(
const model &md,
size_type,
5576 const model::varnamelist &vl,
5577 const model::varnamelist &dl,
5578 const model::mimlist &mims,
5579 model::complex_matlist &matl,
5580 model::complex_veclist &,
5581 model::complex_veclist &,
5583 build_version)
const {
5584 GMM_ASSERT1(matl.size() == 1,
5585 "Fourier-Robin brick has one and only one term");
5586 GMM_ASSERT1(mims.size() == 1,
5587 "Fourier-Robin brick need one and only one mesh_im");
5588 GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
5589 "Wrong number of variables for Fourier-Robin brick");
5591 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5592 const mesh &m = mf_u.linked_mesh();
5594 const mesh_im &mim = *mims[0];
5595 const mesh_fem *mf_a = 0;
5596 mesh_region rg(region);
5597 m.intersect_with_mpi_region(rg);
5598 const model_complex_plain_vector *A = &(md.complex_variable(dl[0]));
5599 mf_a = md.pmesh_fem_of_variable(dl[0]);
5600 s = gmm::vect_size(*A);
5601 if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
5602 GMM_ASSERT1(s == Q*Q,
5603 "Bad format Fourier-Robin brick coefficient");
5605 GMM_TRACE2(
"Fourier-Robin term assembly");
5610 asm_homogeneous_qu_term(matl[0], mim, mf_u, *A, rg);
5613 Fourier_Robin_brick() {
5614 set_flags(
"Fourier Robin condition",
true ,
5623 const std::string &varname,
5624 const std::string &dataexpr,
5627 pbrick pbr = std::make_shared<Fourier_Robin_brick>();
5629 tl.push_back(model::term_description(varname, varname,
true));
5630 return md.
add_brick(pbr, model::varnamelist(1, varname),
5631 model::varnamelist(1, dataexpr), tl,
5632 model::mimlist(1, &mim), region);
5634 std::string test_varname
5635 =
"Test_" + sup_previous_and_dot_to_varname(varname);
5636 std::string expr =
"(("+dataexpr+
")*"+varname+
")."+test_varname;
5638 "Fourier-Robin",
true);
5641 "Fourier-Robin (nonlinear)");
5652 struct have_private_data_brick :
public virtual_brick {
5654 model_real_sparse_matrix rB;
5655 model_complex_sparse_matrix cB;
5656 model_real_plain_vector rL;
5657 model_complex_plain_vector cL;
5661 struct constraint_brick :
public have_private_data_brick {
5663 virtual void real_pre_assembly_in_serial(
const model &md,
size_type,
5664 const model::varnamelist &vl,
5665 const model::varnamelist &dl,
5666 const model::mimlist &mims,
5667 model::real_matlist &matl,
5668 model::real_veclist &vecl,
5669 model::real_veclist &,
5671 if (MPI_IS_MASTER()) {
5673 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5674 "Constraint brick has one and only one term");
5675 GMM_ASSERT1(mims.size() == 0,
5676 "Constraint brick need no mesh_im");
5677 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() <= 1,
5678 "Wrong number of variables for constraint brick");
5680 bool penalized = (vl.size() == 1);
5681 const model_real_plain_vector *COEFF = 0;
5683 bool has_data = (nameL.compare(
"") != 0);
5685 GMM_ASSERT1(nameL.compare(dl.back()) == 0 &&
5686 md.variable_exists(nameL) && md.is_data(nameL),
5688 const model_real_plain_vector &
5689 rrL = has_data ? md.real_variable(nameL) : rL;
5692 COEFF = &(md.real_variable(dl[0]));
5693 GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
5694 "Data for coefficient should be a scalar");
5696 gmm::mult(gmm::transposed(rB),
5697 gmm::scaled(rrL, gmm::abs((*COEFF)[0])), vecl[0]);
5698 gmm::mult(gmm::transposed(rB),
5699 gmm::scaled(rB, gmm::abs((*COEFF)[0])), matl[0]);
5701 gmm::copy(rrL, vecl[0]);
5702 gmm::copy(rB, matl[0]);
5707 virtual void complex_pre_assembly_in_serial(
const model &md,
size_type,
5708 const model::varnamelist &vl,
5709 const model::varnamelist &dl,
5710 const model::mimlist &mims,
5711 model::complex_matlist &matl,
5712 model::complex_veclist &vecl,
5713 model::complex_veclist &,
5715 build_version)
const {
5716 if (MPI_IS_MASTER()) {
5718 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5719 "Constraint brick has one and only one term");
5720 GMM_ASSERT1(mims.size() == 0,
5721 "Constraint brick need no mesh_im");
5722 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() <= 1,
5723 "Wrong number of variables for constraint brick");
5725 bool penalized = (vl.size() == 1);
5726 const model_complex_plain_vector *COEFF = 0;
5728 bool has_data = (nameL.compare(
"") != 0);
5730 GMM_ASSERT1(nameL.compare(dl.back()) == 0 &&
5731 md.variable_exists(nameL) && md.is_data(nameL),
5733 const model_complex_plain_vector &
5734 ccL = has_data ? md.complex_variable(nameL) : cL;
5737 COEFF = &(md.complex_variable(dl[0]));
5738 GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
5739 "Data for coefficient should be a scalar");
5741 gmm::mult(gmm::transposed(cB),
5742 gmm::scaled(ccL, gmm::abs((*COEFF)[0])), vecl[0]);
5743 gmm::mult(gmm::transposed(cB),
5744 gmm::scaled(cB, gmm::abs((*COEFF)[0])), matl[0]);
5746 gmm::copy(ccL, vecl[0]);
5747 gmm::copy(cB, matl[0]);
5752 virtual std::string declare_volume_assembly_string
5753 (
const model &,
size_type,
const model::varnamelist &,
5754 const model::varnamelist &)
const {
5755 return std::string();
5758 constraint_brick(
bool penalized) {
5759 set_flags(penalized ?
"Constraint with penalization brick"
5760 :
"Constraint with multipliers brick",
5769 model_real_sparse_matrix &set_private_data_brick_real_matrix
5771 pbrick pbr = md.brick_pointer(indbrick);
5772 md.touch_brick(indbrick);
5773 have_private_data_brick *p =
dynamic_cast<have_private_data_brick *
>
5774 (
const_cast<virtual_brick *
>(pbr.get()));
5775 GMM_ASSERT1(p,
"Wrong type of brick");
5779 model_real_plain_vector &set_private_data_brick_real_rhs
5781 pbrick pbr = md.brick_pointer(indbrick);
5782 md.touch_brick(indbrick);
5783 have_private_data_brick *p =
dynamic_cast<have_private_data_brick *
>
5784 (
const_cast<virtual_brick *
>(pbr.get()));
5785 GMM_ASSERT1(p,
"Wrong type of brick");
5786 if (p->nameL.compare(
"") != 0) GMM_WARNING1(
"Rhs already set by data name");
5790 model_complex_sparse_matrix &set_private_data_brick_complex_matrix
5792 pbrick pbr = md.brick_pointer(indbrick);
5793 md.touch_brick(indbrick);
5794 have_private_data_brick *p =
dynamic_cast<have_private_data_brick *
>
5795 (
const_cast<virtual_brick *
>(pbr.get()));
5796 GMM_ASSERT1(p,
"Wrong type of brick");
5800 model_complex_plain_vector &set_private_data_brick_complex_rhs
5802 pbrick pbr = md.brick_pointer(indbrick);
5803 md.touch_brick(indbrick);
5804 have_private_data_brick *p =
dynamic_cast<have_private_data_brick *
>
5805 (
const_cast<virtual_brick *
>(pbr.get()));
5806 GMM_ASSERT1(p,
"Wrong type of brick");
5807 if (p->nameL.compare(
"") != 0) GMM_WARNING1(
"Rhs already set by data name");
5811 void set_private_data_rhs
5812 (model &md,
size_type indbrick,
const std::string &varname) {
5813 pbrick pbr = md.brick_pointer(indbrick);
5814 md.touch_brick(indbrick);
5815 have_private_data_brick *p =
dynamic_cast<have_private_data_brick *
>
5816 (
const_cast<virtual_brick *
>(pbr.get()));
5817 GMM_ASSERT1(p,
"Wrong type of brick");
5818 if (p->nameL.compare(varname) != 0) {
5819 model::varnamelist dl = md.datanamelist_of_brick(indbrick);
5820 if (p->nameL.compare(
"") == 0) dl.push_back(varname);
5821 else dl.back() = varname;
5822 md.change_data_of_brick(indbrick, dl);
5827 size_type add_constraint_with_penalization
5828 (model &md,
const std::string &varname, scalar_type penalisation_coeff) {
5829 std::string coeffname = md.new_name(
"penalization_on_" + varname);
5830 md.add_fixed_size_data(coeffname, 1);
5831 if (md.is_complex())
5832 md.set_complex_variable(coeffname)[0] = penalisation_coeff;
5834 md.set_real_variable(coeffname)[0] = penalisation_coeff;
5835 pbrick pbr = std::make_shared<constraint_brick>(
true);
5837 tl.push_back(model::term_description(varname, varname,
true));
5838 model::varnamelist vl(1, varname);
5839 model::varnamelist dl(1, coeffname);
5840 return md.add_brick(pbr, vl, dl, tl, model::mimlist(),
size_type(-1));
5843 size_type add_constraint_with_multipliers
5844 (model &md,
const std::string &varname,
const std::string &multname) {
5845 pbrick pbr = std::make_shared<constraint_brick>(
false);
5847 tl.push_back(model::term_description(multname, varname,
true));
5848 model::varnamelist vl(1, varname);
5849 vl.push_back(multname);
5850 model::varnamelist dl;
5851 return md.add_brick(pbr, vl, dl, tl, model::mimlist(),
size_type(-1));
5861 struct explicit_matrix_brick :
public have_private_data_brick {
5863 virtual void real_pre_assembly_in_serial(
const model &,
size_type,
5864 const model::varnamelist &vl,
5865 const model::varnamelist &dl,
5866 const model::mimlist &mims,
5867 model::real_matlist &matl,
5868 model::real_veclist &vecl,
5869 model::real_veclist &,
5871 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5872 "Explicit matrix has one and only one term");
5873 GMM_ASSERT1(mims.size() == 0,
"Explicit matrix need no mesh_im");
5874 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() == 0,
5875 "Wrong number of variables for explicit matrix brick");
5876 GMM_ASSERT1(gmm::mat_ncols(rB) == gmm::mat_ncols(matl[0]) &&
5877 gmm::mat_nrows(rB) == gmm::mat_nrows(matl[0]),
5878 "Explicit matrix brick dimension mismatch ("<<
5879 gmm::mat_ncols(rB)<<
"x"<<gmm::mat_nrows(rB)<<
") != ("<<
5880 gmm::mat_ncols(matl[0])<<
"x"<<gmm::mat_nrows(matl[0])<<
")");
5881 gmm::copy(rB, matl[0]);
5884 virtual void complex_pre_assembly_in_serial(
const model &,
size_type,
5885 const model::varnamelist &vl,
5886 const model::varnamelist &dl,
5887 const model::mimlist &mims,
5888 model::complex_matlist &matl,
5889 model::complex_veclist &vecl,
5890 model::complex_veclist &,
5892 build_version)
const {
5893 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5894 "Explicit matrix has one and only one term");
5895 GMM_ASSERT1(mims.size() == 0,
"Explicit matrix need no mesh_im");
5896 GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() == 0,
5897 "Wrong number of variables for explicit matrix brick");
5898 gmm::copy(cB, matl[0]);
5901 virtual std::string declare_volume_assembly_string
5902 (
const model &,
size_type,
const model::varnamelist &,
5903 const model::varnamelist &)
const {
5904 return std::string();
5907 explicit_matrix_brick(
bool symmetric_,
bool coercive_) {
5908 set_flags(
"Explicit matrix brick",
5910 symmetric_ , coercive_ ,
5917 (model &md,
const std::string &varname1,
const std::string &varname2,
5918 bool issymmetric,
bool iscoercive) {
5919 pbrick pbr = std::make_shared<explicit_matrix_brick>(issymmetric,
5922 tl.push_back(model::term_description(varname1, varname2, issymmetric));
5923 model::varnamelist vl(1, varname1);
5924 vl.push_back(varname2);
5925 model::varnamelist dl;
5926 return md.add_brick(pbr, vl, dl, tl, model::mimlist(),
size_type(-1));
5935 struct explicit_rhs_brick :
public have_private_data_brick {
5937 virtual void real_pre_assembly_in_serial(
const model &,
size_type,
5938 const model::varnamelist &vl,
5939 const model::varnamelist &dl,
5940 const model::mimlist &mims,
5941 model::real_matlist &matl,
5942 model::real_veclist &vecl,
5943 model::real_veclist &,
5945 if (MPI_IS_MASTER()) {
5946 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5947 "Explicit rhs has one and only one term");
5948 GMM_ASSERT1(mims.size() == 0,
"Explicit rhs need no mesh_im");
5949 GMM_ASSERT1(vl.size() == 1 && dl.size() == 0,
5950 "Wrong number of variables for explicit rhs brick");
5951 gmm::copy(rL, vecl[0]);
5955 virtual void complex_pre_assembly_in_serial(
const model &,
size_type,
5956 const model::varnamelist &vl,
5957 const model::varnamelist &dl,
5958 const model::mimlist &mims,
5959 model::complex_matlist &matl,
5960 model::complex_veclist &vecl,
5961 model::complex_veclist &,
5963 build_version)
const {
5964 if (MPI_IS_MASTER()) {
5965 GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5966 "Explicit rhs has one and only one term");
5967 GMM_ASSERT1(mims.size() == 0,
"Explicit rhs need no mesh_im");
5968 GMM_ASSERT1(vl.size() == 1 && dl.size() == 0,
5969 "Wrong number of variables for explicit rhs brick");
5970 gmm::copy(cL, vecl[0]);
5975 virtual std::string declare_volume_assembly_string
5976 (
const model &,
size_type,
const model::varnamelist &,
5977 const model::varnamelist &)
const {
5978 return std::string();
5981 explicit_rhs_brick() {
5982 set_flags(
"Explicit rhs brick",
5992 (model &md,
const std::string &varname) {
5993 pbrick pbr = std::make_shared<explicit_rhs_brick>();
5995 tl.push_back(model::term_description(varname));
5996 model::varnamelist vl(1, varname);
5997 model::varnamelist dl;
5998 return md.add_brick(pbr, vl, dl, tl, model::mimlist(),
size_type(-1));
6008 struct iso_lin_elasticity_new_brick :
public virtual_brick {
6010 std::string expr, dataname3;
6012 void asm_real_tangent_terms(
const model &md,
size_type ib,
6013 const model::varnamelist &vl,
6014 const model::varnamelist &dl,
6015 const model::mimlist &mims,
6016 model::real_matlist &matl,
6017 model::real_veclist &vecl,
6018 model::real_veclist &,
6020 build_version version)
const override {
6021 GMM_ASSERT1(vl.size() == 1,
"Linearized isotropic elasticity brick "
6022 "has one and only one variable");
6023 GMM_ASSERT1(matl.size() == 1,
"Linearized isotropic elasticity brick "
6024 "has one and only one term");
6025 GMM_ASSERT1(mims.size() == 1,
"Linearized isotropic elasticity brick "
6026 "needs one and only one mesh_im");
6028 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0);
6029 for (
size_type i = 0; i < dl.size(); ++i) {
6030 recompute_matrix = recompute_matrix ||
6031 md.is_var_newer_than_brick(dl[i], ib);
6034 if (recompute_matrix) {
6036 ga_workspace workspace(md, ga_workspace::inherit::ALL);
6037 workspace.add_expression(expr, *(mims[0]), region);
6038 GMM_TRACE2(name <<
": generic matrix assembly");
6039 workspace.assembly(2);
6040 scalar_type
alpha = scalar_type(1)
6041 / (workspace.factor_of_variable(vl[0]));
6042 const auto &R=workspace.assembled_matrix();
6043 gmm::sub_interval I = workspace.interval_of_variable(vl[0]);
6044 gmm::copy(gmm::scaled(gmm::sub_matrix(R, I, I), alpha),
6048 if (dataname3.size()) {
6054 gmm::scaled(md.real_variable(dataname3), scalar_type(-1)),
6060 void real_post_assembly_in_serial(
const model &md,
size_type ib,
6061 const model::varnamelist &,
6062 const model::varnamelist &,
6063 const model::mimlist &,
6064 model::real_matlist &,
6065 model::real_veclist &vecl,
6066 model::real_veclist &,
6068 build_version)
const override {
6069 md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
6073 virtual std::string declare_volume_assembly_string
6074 (
const model &,
size_type,
const model::varnamelist &,
6075 const model::varnamelist &)
const {
6079 iso_lin_elasticity_new_brick(
const std::string &expr_,
6080 const std::string &dataname3_) {
6081 expr = expr_; dataname3 = dataname3_;
6082 set_flags(
"Linearized isotropic elasticity",
true ,
6092 const std::string &dataexpr1,
const std::string &dataexpr2,
6093 size_type region,
const std::string &dataname3) {
6094 std::string test_varname
6095 =
"Test_" + sup_previous_and_dot_to_varname(varname);
6097 std::string expr1 =
"((("+dataexpr1+
")*(Div_"+varname+
"-Div_"+dataname3
6098 +
"))*Id(meshdim)+(2*("+dataexpr2+
"))*(Sym(Grad_"+varname
6099 +
")-Sym(Grad_"+dataname3+
"))):Grad_" +test_varname;
6100 std::string expr2 =
"(Div_"+varname+
"*(("+dataexpr1+
")*Id(meshdim))"
6101 +
"+(2*("+dataexpr2+
"))*Sym(Grad_"+varname+
")):Grad_"+test_varname;
6104 model::varnamelist vl, dl;
6106 ga_workspace workspace(md, ga_workspace::inherit::ALL);
6107 workspace.add_expression(expr2, mim, region);
6108 model::varnamelist vl_test1, vl_test2;
6109 is_lin = workspace.used_variables(vl, vl_test1, vl_test2, dl, 2);
6112 pbrick pbr = std::make_shared<iso_lin_elasticity_new_brick>
6115 tl.push_back(model::term_description(varname,
6116 sup_previous_and_dot_to_varname(varname),
true));
6117 if (dataname3.size()) dl.push_back(dataname3);
6118 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
6120 return add_nonlinear_generic_assembly_brick
6121 (md, mim, dataname3.size() ? expr1 : expr2, region,
false,
false,
6122 "Linearized isotropic elasticity (with nonlinear dependance)");
6128 const std::string &data_E,
const std::string &data_nu,
6130 std::string test_varname
6131 =
"Test_" + sup_previous_and_dot_to_varname(varname);
6133 std::string mu =
"(("+data_E+
")/(2*(1+("+data_nu+
"))))";
6134 std::string lambda =
"(("+data_E+
")*("+data_nu+
")/((1+("+data_nu+
"))*(1-2*("
6136 std::string expr = lambda+
"*Div_"+varname+
"*Div_"+test_varname
6137 +
"+"+mu+
"*(Grad_"+varname+
"+Grad_"+varname+
"'):Grad_"+test_varname;
6141 ga_workspace workspace(md, ga_workspace::inherit::ALL);
6142 workspace.add_expression(expr, mim, region);
6143 is_lin = workspace.is_linear(2);
6147 "Linearized isotropic elasticity");
6150 (md, mim, expr, region,
false,
false,
6151 "Linearized isotropic elasticity (with nonlinear dependance)");
6157 const std::string &data_E,
const std::string &data_nu,
6159 std::string test_varname
6160 =
"Test_" + sup_previous_and_dot_to_varname(varname);
6163 GMM_ASSERT1(mfu,
"The variable should be a fem variable");
6166 std::string mu =
"(("+data_E+
")/(2*(1+("+data_nu+
"))))";
6167 std::string lambda =
"(("+data_E+
")*("+data_nu+
")/((1+("+data_nu+
"))*(1-2*("
6170 lambda =
"(("+data_E+
")*("+data_nu+
")/((1-sqr("+data_nu+
"))))";
6171 std::string expr = lambda+
"*Div_"+varname+
"*Div_"+test_varname
6172 +
"+"+mu+
"*(Grad_"+varname+
"+Grad_"+varname+
"'):Grad_"+test_varname;
6176 ga_workspace workspace(md, ga_workspace::inherit::ALL);
6177 workspace.add_expression(expr, mim, region);
6178 is_lin = workspace.is_linear(2);
6182 "Linearized isotropic elasticity");
6185 (md, mim, expr, region,
false,
false,
6186 "Linearized isotropic elasticity (with nonlinear dependance)");
6191 void compute_isotropic_linearized_Von_Mises_or_Tresca
6192 (model &md,
const std::string &varname,
const std::string &data_lambda,
6193 const std::string &data_mu,
const mesh_fem &mf_vm,
6194 model_real_plain_vector &VM,
bool tresca) {
6197 const mesh_fem &mf_u = md.mesh_fem_of_variable(varname);
6198 const mesh_fem *mf_lambda = md.pmesh_fem_of_variable(data_lambda);
6199 const model_real_plain_vector *lambda=&(md.real_variable(data_lambda));
6200 const mesh_fem *mf_mu = md.pmesh_fem_of_variable(data_mu);
6201 const model_real_plain_vector *mu = &(md.real_variable(data_mu));
6204 if (mf_lambda) sl = sl * mf_lambda->get_qdim() / mf_lambda->nb_dof();
6206 if (mf_mu) sm = sm * mf_mu->get_qdim() / mf_mu->nb_dof();
6208 GMM_ASSERT1(sl == 1 && sm == 1,
"Bad format for Lame coefficients");
6209 GMM_ASSERT1(mf_lambda == mf_mu,
6210 "The two Lame coefficients should be described on the same "
6211 "finite element method.");
6215 md.real_variable(varname), VM,
6216 *mf_lambda, *lambda,
6221 model_real_plain_vector LAMBDA(mf_lambda->nb_dof(), (*lambda)[0]);
6222 model_real_plain_vector MU(mf_lambda->nb_dof(), (*mu)[0]);
6224 md.real_variable(varname), VM,
6233 std::string sigma_d=
"("+data_mu+
")*(Grad_"+varname+
"+Grad_"+varname+
"')";
6234 std::string expr =
"sqrt(3/2)*Norm(Deviator("+sigma_d+
"))";
6235 ga_interpolation_Lagrange_fem(md, expr, mf_vm, VM);
6241 (
model &md,
const std::string &varname,
const std::string &data_E,
6242 const std::string &data_nu,
const mesh_fem &mf_vm,
6243 model_real_plain_vector &VM) {
6247 std::string mu =
"(("+data_E+
")/(2*(1+("+data_nu+
"))))";
6250 std::string sigma_d = mu+
"*(Grad_"+varname+
"+Grad_"+varname+
"')";
6251 std::string expr =
"sqrt(3/2)*Norm(Deviator("+sigma_d+
"))";
6252 ga_interpolation_Lagrange_fem(md, expr, mf_vm, VM);
6256 (
model &md,
const std::string &varname,
const std::string &data_E,
6257 const std::string &data_nu,
const mesh_fem &mf_vm,
6258 model_real_plain_vector &VM) {
6267 std::string mu =
"(("+data_E+
")/(2*(1+("+data_nu+
"))))";
6270 std::string sigma_d = mu+
"*(Grad_"+varname+
"+Grad_"+varname+
"')";
6271 std::string expr =
"sqrt(3/2)*Norm(Deviator("+sigma_d+
"))";
6272 ga_interpolation_Lagrange_fem(md, expr, mf_vm, VM);
6282 struct linear_incompressibility_brick :
public virtual_brick {
6284 virtual void asm_real_tangent_terms(
const model &md,
size_type ,
6285 const model::varnamelist &vl,
6286 const model::varnamelist &dl,
6287 const model::mimlist &mims,
6288 model::real_matlist &matl,
6289 model::real_veclist &,
6290 model::real_veclist &,
6292 build_version)
const {
6294 GMM_ASSERT1((matl.size() == 1 && dl.size() == 0)
6295 || (matl.size() == 2 && dl.size() == 1),
6296 "Wrong term and/or data number for Linear incompressibility "
6298 GMM_ASSERT1(mims.size() == 1,
"Linear incompressibility brick need one "
6299 "and only one mesh_im");
6300 GMM_ASSERT1(vl.size() == 2,
"Wrong number of variables for linear "
6301 "incompressibility brick");
6303 bool penalized = (dl.size() == 1);
6304 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6305 const mesh_fem &mf_p = md.mesh_fem_of_variable(vl[1]);
6306 const mesh_im &mim = *mims[0];
6307 const model_real_plain_vector *COEFF = 0;
6308 const mesh_fem *mf_data = 0;
6311 COEFF = &(md.real_variable(dl[0]));
6312 mf_data = md.pmesh_fem_of_variable(dl[0]);
6314 if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
6315 GMM_ASSERT1(s == 1,
"Bad format for the penalization parameter");
6318 mesh_region rg(region);
6319 mim.linked_mesh().intersect_with_mpi_region(rg);
6321 GMM_TRACE2(
"Stokes term assembly");
6329 gmm::scale(matl[1], scalar_type(-1));
6333 gmm::scale(matl[1], -(*COEFF)[0]);
6340 virtual void real_post_assembly_in_serial(
const model &,
size_type,
6341 const model::varnamelist &,
6342 const model::varnamelist &,
6343 const model::mimlist &,
6344 model::real_matlist &,
6345 model::real_veclist &,
6346 model::real_veclist &,
6348 build_version)
const
6352 linear_incompressibility_brick() {
6353 set_flags(
"Linear incompressibility brick",
6363 const std::string &multname,
size_type region,
6364 const std::string &dataexpr) {
6366 pbrick pbr = std::make_shared<linear_incompressibility_brick>();
6368 tl.push_back(model::term_description(multname, varname,
true));
6369 model::varnamelist vl(1, varname);
6370 vl.push_back(multname);
6371 model::varnamelist dl;
6372 if (dataname.size()) {
6373 dl.push_back(dataexpr);
6374 tl.push_back(model::term_description(multname, multname,
true));
6376 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
6378 std::string test_varname
6379 =
"Test_" + sup_previous_and_dot_to_varname(varname);
6380 std::string test_multname
6381 =
"Test_" + sup_previous_and_dot_to_varname(multname);
6383 if (dataexpr.size())
6384 expr =
"-"+multname+
"*Div_"+test_varname +
"-"+test_multname
6385 +
"*Div_"+varname+
"+(("+dataexpr+
")*"+multname+
")*"+test_multname;
6387 expr =
"-"+multname+
"*Div_"+test_varname +
"-"+test_multname
6390 "Linear incompressibility",
true);
6393 (md, mim, expr, region,
false,
false,
6394 "Linear incompressibility (with nonlinear dependance)");
6407 struct mass_brick :
public virtual_brick {
6409 virtual void asm_real_tangent_terms(
const model &md,
size_type,
6410 const model::varnamelist &vl,
6411 const model::varnamelist &dl,
6412 const model::mimlist &mims,
6413 model::real_matlist &matl,
6414 model::real_veclist &,
6415 model::real_veclist &,
6417 build_version)
const {
6418 GMM_ASSERT1(matl.size() == 1,
6419 "Mass brick has one and only one term");
6420 GMM_ASSERT1(mims.size() == 1,
6421 "Mass brick need one and only one mesh_im");
6422 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
6423 "Wrong number of variables for mass brick");
6425 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6426 const mesh &m = mf_u.linked_mesh();
6427 const mesh_im &mim = *mims[0];
6428 mesh_region rg(region);
6429 m.intersect_with_mpi_region(rg);
6431 const mesh_fem *mf_rho = 0;
6432 const model_real_plain_vector *rho = 0;
6435 mf_rho = md.pmesh_fem_of_variable(dl[0]);
6436 rho = &(md.real_variable(dl[0]));
6438 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6439 GMM_ASSERT1(sl == 1,
"Bad format of mass brick coefficient");
6442 GMM_TRACE2(
"Mass matrix assembly");
6444 if (dl.size() && mf_rho) {
6448 if (dl.size()) gmm::scale(matl[0], (*rho)[0]);
6452 virtual void asm_complex_tangent_terms(
const model &md,
size_type,
6453 const model::varnamelist &vl,
6454 const model::varnamelist &dl,
6455 const model::mimlist &mims,
6456 model::complex_matlist &matl,
6457 model::complex_veclist &,
6458 model::complex_veclist &,
6460 build_version)
const {
6461 GMM_ASSERT1(matl.size() == 1,
6462 "Mass brick has one and only one term");
6463 GMM_ASSERT1(mims.size() == 1,
6464 "Mass brick need one and only one mesh_im");
6465 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
6466 "Wrong number of variables for mass brick");
6468 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6469 const mesh &m = mf_u.linked_mesh();
6470 const mesh_im &mim = *mims[0];
6471 mesh_region rg(region);
6472 m.intersect_with_mpi_region(rg);
6474 const mesh_fem *mf_rho = 0;
6475 const model_complex_plain_vector *rho = 0;
6478 mf_rho = md.pmesh_fem_of_variable(dl[0]);
6479 rho = &(md.complex_variable(dl[0]));
6481 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6482 GMM_ASSERT1(sl == 1,
"Bad format of mass brick coefficient");
6485 GMM_TRACE2(
"Mass matrix assembly");
6487 if (dl.size() && mf_rho) {
6491 if (dl.size()) gmm::scale(matl[0], (*rho)[0]);
6495 virtual std::string declare_volume_assembly_string
6496 (
const model &,
size_type,
const model::varnamelist &,
6497 const model::varnamelist &)
const {
6498 return std::string();
6502 set_flags(
"Mass brick",
true ,
6512 const std::string &dataexpr_rho,
size_type region) {
6514 pbrick pbr = std::make_shared<mass_brick>();
6516 tl.push_back(model::term_description(varname, varname,
true));
6517 model::varnamelist dl;
6518 if (dataexpr_rho.size())
6519 dl.push_back(dataexpr_rho);
6520 return md.
add_brick(pbr, model::varnamelist(1, varname), dl, tl,
6521 model::mimlist(1, &mim), region);
6523 std::string test_varname
6524 =
"Test_" + sup_previous_and_dot_to_varname(varname);
6526 if (dataexpr_rho.size())
6527 expr =
"(("+dataexpr_rho+
")*"+varname+
")."+test_varname;
6529 expr = varname+
"."+test_varname;
6531 "Mass matrix",
true);
6534 "Mass matrix (nonlinear)");
6545 struct lumped_mass_brick_for_first_order :
public virtual_brick {
6547 virtual void asm_real_tangent_terms(
const model &md,
size_type,
6548 const model::varnamelist &vl,
6549 const model::varnamelist &dl,
6550 const model::mimlist &mims,
6551 model::real_matlist &matl,
6552 model::real_veclist &,
6553 model::real_veclist &,
6555 build_version)
const {
6556 GMM_ASSERT1(matl.size() == 1,
6557 "Lumped Mass brick has one and only one term");
6558 GMM_ASSERT1(mims.size() == 1,
6559 "Lumped Mass brick needs one and only one mesh_im");
6560 GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
6561 "Wrong number of variables for lumped mass brick");
6563 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6564 const mesh &m = mf_u.linked_mesh();
6565 const mesh_im &mim = *mims[0];
6566 mesh_region rg(region);
6567 m.intersect_with_mpi_region(rg);
6569 const mesh_fem *mf_rho = 0;
6570 const model_real_plain_vector *rho = 0;
6573 mf_rho = md.pmesh_fem_of_variable(dl[0]);
6574 rho = &(md.real_variable(dl[0]));
6576 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6577 GMM_ASSERT1(sl == 1,
"Bad format of mass brick coefficient");
6580 GMM_TRACE2(
"Lumped mass matrix assembly (please check that integration is 1st order.)");
6582 if (dl.size() && mf_rho) {
6586 if (dl.size()) gmm::scale(matl[0], (*rho)[0]);
6591 lumped_mass_brick_for_first_order() {
6592 set_flags(
"Lumped mass brick",
true ,
6602 const std::string &dataexpr_rho,
size_type region) {
6603 pbrick pbr = std::make_shared<lumped_mass_brick_for_first_order>();
6605 tl.push_back(model::term_description(varname, varname,
true));
6606 model::varnamelist dl;
6607 if (dataexpr_rho.size())
6608 dl.push_back(dataexpr_rho);
6609 return md.
add_brick(pbr, model::varnamelist(1, varname), dl, tl,
6610 model::mimlist(1, &mim), region);
6626 struct basic_d_on_dt_brick :
public virtual_brick {
6628 virtual void asm_real_tangent_terms(
const model &md,
size_type ib,
6629 const model::varnamelist &vl,
6630 const model::varnamelist &dl,
6631 const model::mimlist &mims,
6632 model::real_matlist &matl,
6633 model::real_veclist &vecl,
6634 model::real_veclist &,
6636 build_version version)
const {
6637 GMM_ASSERT1(matl.size() == 1,
6638 "Basic d/dt brick has one and only one term");
6639 GMM_ASSERT1(mims.size() == 1,
6640 "Basic d/dt brick need one and only one mesh_im");
6641 GMM_ASSERT1(vl.size() == 1 && dl.size() >= 2 && dl.size() <= 3,
6642 "Wrong number of variables for basic d/dt brick");
6646 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
6647 || (md.is_var_newer_than_brick(dl[1], ib));
6649 recompute_matrix = recompute_matrix ||
6650 md.is_var_newer_than_brick(dl[2], ib);
6653 if (recompute_matrix) {
6654 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6655 const mesh &m = mf_u.linked_mesh();
6656 const mesh_im &mim = *mims[0];
6657 mesh_region rg(region);
6658 m.intersect_with_mpi_region(rg);
6660 const model_real_plain_vector &dt = md.real_variable(dl[1]);
6661 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for time step");
6663 const mesh_fem *mf_rho = 0;
6664 const model_real_plain_vector *rho = 0;
6666 if (dl.size() > 2) {
6667 mf_rho = md.pmesh_fem_of_variable(dl[2]);
6668 rho = &(md.real_variable(dl[2]));
6670 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6671 GMM_ASSERT1(sl == 1,
"Bad format for density");
6674 GMM_TRACE2(
"Mass matrix assembly for d_on_dt brick");
6675 if (dl.size() > 2 && mf_rho) {
6678 gmm::scale(matl[0], scalar_type(1) / dt[0]);
6682 if (dl.size() > 2) gmm::scale(matl[0], (*rho)[0] / dt[0]);
6683 else gmm::scale(matl[0], scalar_type(1) / dt[0]);
6686 gmm::mult(matl[0], md.real_variable(dl[0], 1), vecl[0]);
6689 virtual void asm_complex_tangent_terms(
const model &md,
size_type ib,
6690 const model::varnamelist &vl,
6691 const model::varnamelist &dl,
6692 const model::mimlist &mims,
6693 model::complex_matlist &matl,
6694 model::complex_veclist &vecl,
6695 model::complex_veclist &,
6697 build_version version)
const {
6698 GMM_ASSERT1(matl.size() == 1,
6699 "Basic d/dt brick has one and only one term");
6700 GMM_ASSERT1(mims.size() == 1,
6701 "Basic d/dt brick need one and only one mesh_im");
6702 GMM_ASSERT1(vl.size() == 1 && dl.size() >= 2 && dl.size() <= 3,
6703 "Wrong number of variables for basic d/dt brick");
6707 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
6708 || (md.is_var_newer_than_brick(dl[1], ib));
6710 recompute_matrix = recompute_matrix ||
6711 md.is_var_newer_than_brick(dl[2], ib);
6713 if (recompute_matrix) {
6714 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6715 const mesh &m = mf_u.linked_mesh();
6716 const mesh_im &mim = *mims[0];
6718 mesh_region rg(region);
6719 m.intersect_with_mpi_region(rg);
6721 const model_complex_plain_vector &dt = md.complex_variable(dl[1]);
6722 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for time step");
6724 const mesh_fem *mf_rho = 0;
6725 const model_complex_plain_vector *rho = 0;
6727 if (dl.size() > 2) {
6728 mf_rho = md.pmesh_fem_of_variable(dl[2]);
6729 rho = &(md.complex_variable(dl[2]));
6731 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6732 GMM_ASSERT1(sl == 1,
"Bad format for density");
6735 GMM_TRACE2(
"Mass matrix assembly for d_on_dt brick");
6736 if (dl.size() > 2 && mf_rho) {
6739 gmm::scale(matl[0], scalar_type(1) / dt[0]);
6743 if (dl.size() > 2) gmm::scale(matl[0], (*rho)[0] / dt[0]);
6744 else gmm::scale(matl[0], scalar_type(1) / dt[0]);
6747 gmm::mult(matl[0], md.complex_variable(dl[0], 1), vecl[0]);
6750 virtual std::string declare_volume_assembly_string
6751 (
const model &,
size_type,
const model::varnamelist &,
6752 const model::varnamelist &)
const {
6753 return std::string();
6756 basic_d_on_dt_brick() {
6757 set_flags(
"Basic d/dt brick",
true ,
6767 const std::string &dataname_dt,
const std::string &dataname_rho,
6769 pbrick pbr = std::make_shared<basic_d_on_dt_brick>();
6771 tl.push_back(model::term_description(varname, varname,
true));
6772 model::varnamelist dl(1, varname);
6773 dl.push_back(dataname_dt);
6774 if (dataname_rho.size())
6775 dl.push_back(dataname_rho);
6776 return md.
add_brick(pbr, model::varnamelist(1, varname), dl, tl,
6777 model::mimlist(1, &mim), region);
6788 struct basic_d2_on_dt2_brick :
public virtual_brick {
6790 mutable scalar_type old_alphadt2;
6792 virtual void asm_real_tangent_terms(
const model &md,
size_type ib,
6793 const model::varnamelist &vl,
6794 const model::varnamelist &dl,
6795 const model::mimlist &mims,
6796 model::real_matlist &matl,
6797 model::real_veclist &vecl,
6798 model::real_veclist &,
6800 build_version version)
const {
6801 GMM_ASSERT1(matl.size() == 1,
6802 "Basic d2/dt2 brick has one and only one term");
6803 GMM_ASSERT1(mims.size() == 1,
6804 "Basic d2/dt2 brick need one and only one mesh_im");
6805 GMM_ASSERT1(vl.size() == 1 && dl.size() >= 4 && dl.size() <= 5,
6806 "Wrong number of variables for basic d2/dt2 brick");
6808 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0);
6810 recompute_matrix = recompute_matrix || md.is_var_newer_than_brick(dl[2], ib);
6811 if (dl.size() > 4) recompute_matrix || md.is_var_newer_than_brick(dl[4], ib);
6813 const model_real_plain_vector &dt = md.real_variable(dl[2]);
6814 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for time step");
6815 const model_real_plain_vector &alpha = md.real_variable(dl[3]);
6816 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for parameter alpha");
6817 scalar_type alphadt2 = gmm::sqr(dt[0]) * alpha[0];
6819 if (!recompute_matrix && alphadt2 != old_alphadt2)
6820 gmm::scale(matl[0], old_alphadt2/alphadt2);
6821 old_alphadt2 = alphadt2;
6823 if (recompute_matrix) {
6824 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6825 const mesh &m = mf_u.linked_mesh();
6826 const mesh_im &mim = *mims[0];
6827 mesh_region rg(region);
6828 m.intersect_with_mpi_region(rg);
6830 const mesh_fem *mf_rho = 0;
6831 const model_real_plain_vector *rho = 0;
6833 if (dl.size() > 4) {
6834 mf_rho = md.pmesh_fem_of_variable(dl[4]);
6835 rho = &(md.real_variable(dl[4]));
6837 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6838 GMM_ASSERT1(sl == 1,
"Bad format for density");
6841 GMM_TRACE2(
"Mass matrix assembly for d2_on_dt2 brick");
6842 if (dl.size() > 4 && mf_rho) {
6845 gmm::scale(matl[0], scalar_type(1) / alphadt2);
6849 if (dl.size() > 4) gmm::scale(matl[0], (*rho)[0] / alphadt2);
6850 else gmm::scale(matl[0], scalar_type(1) / alphadt2);
6853 gmm::mult(matl[0], md.real_variable(dl[0], 1), vecl[0]);
6854 gmm::mult_add(matl[0], gmm::scaled(md.real_variable(dl[1], 1), dt[0]),
6858 virtual void asm_complex_tangent_terms(
const model &md,
size_type ib,
6859 const model::varnamelist &vl,
6860 const model::varnamelist &dl,
6861 const model::mimlist &mims,
6862 model::complex_matlist &matl,
6863 model::complex_veclist &vecl,
6864 model::complex_veclist &,
6866 build_version version)
const {
6867 GMM_ASSERT1(matl.size() == 1,
6868 "Basic d2/dt2 brick has one and only one term");
6869 GMM_ASSERT1(mims.size() == 1,
6870 "Basic d2/dt2 brick need one and only one mesh_im");
6871 GMM_ASSERT1(vl.size() == 1 && dl.size() >= 4 && dl.size() <= 5,
6872 "Wrong number of variables for basic d2/dt2 brick");
6874 bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0);
6876 recompute_matrix = recompute_matrix || md.is_var_newer_than_brick(dl[2], ib);
6877 if (dl.size() > 4) recompute_matrix || md.is_var_newer_than_brick(dl[4], ib);
6880 const model_complex_plain_vector &dt = md.complex_variable(dl[2]);
6881 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for time step");
6882 const model_complex_plain_vector &
alpha = md.complex_variable(dl[3]);
6883 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for parameter alpha");
6884 scalar_type alphadt2 = gmm::real(gmm::sqr(dt[0]) * alpha[0]);
6886 if (!recompute_matrix && alphadt2 != old_alphadt2)
6887 gmm::scale(matl[0], old_alphadt2/alphadt2);
6888 old_alphadt2 = alphadt2;
6890 if (recompute_matrix) {
6891 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6892 const mesh &m = mf_u.linked_mesh();
6893 const mesh_im &mim = *mims[0];
6894 mesh_region rg(region);
6895 m.intersect_with_mpi_region(rg);
6897 const mesh_fem *mf_rho = 0;
6898 const model_complex_plain_vector *rho = 0;
6900 if (dl.size() > 4) {
6901 mf_rho = md.pmesh_fem_of_variable(dl[4]);
6902 rho = &(md.complex_variable(dl[4]));
6904 if (mf_rho) sl = sl * mf_rho->get_qdim() / mf_rho->nb_dof();
6905 GMM_ASSERT1(sl == 1,
"Bad format for density");
6908 GMM_TRACE2(
"Mass matrix assembly for d2_on_dt2 brick");
6909 if (dl.size() > 4 && mf_rho) {
6912 gmm::scale(matl[0], scalar_type(1) / alphadt2);
6916 if (dl.size() > 4) gmm::scale(matl[0], (*rho)[0] / alphadt2);
6917 else gmm::scale(matl[0], scalar_type(1) / alphadt2);
6920 gmm::mult(matl[0], md.complex_variable(dl[0], 1), vecl[0]);
6921 gmm::mult_add(matl[0], gmm::scaled(md.complex_variable(dl[1], 1), dt[0]),
6925 virtual std::string declare_volume_assembly_string
6926 (
const model &,
size_type,
const model::varnamelist &,
6927 const model::varnamelist &)
const {
6928 return std::string();
6931 basic_d2_on_dt2_brick() {
6932 set_flags(
"Basic d2/dt2 brick",
true ,
6942 const std::string &datanameV,
6943 const std::string &dataname_dt,
6944 const std::string &dataname_alpha,
6945 const std::string &dataname_rho,
6947 pbrick pbr = std::make_shared<basic_d2_on_dt2_brick>();
6949 tl.push_back(model::term_description(varnameU, varnameU,
true));
6950 model::varnamelist dl(1, varnameU);
6951 dl.push_back(datanameV);
6952 dl.push_back(dataname_dt);
6953 dl.push_back(dataname_alpha);
6954 if (dataname_rho.size())
6955 dl.push_back(dataname_rho);
6956 return md.
add_brick(pbr, model::varnamelist(1, varnameU), dl, tl,
6957 model::mimlist(1, &mim), region);
6976 void theta_method_dispatcher::set_dispatch_coeff(
const model &md,
size_type ib)
const {
6978 if (md.is_complex())
6979 theta = gmm::real(md.complex_variable(param_names[0])[0]);
6981 theta = md.real_variable(param_names[0])[0];
6983 md.matrix_coeff_of_brick(ib) = theta;
6985 md.rhs_coeffs_of_brick(ib)[0] = theta;
6987 md.rhs_coeffs_of_brick(ib)[1] = (scalar_type(1) - theta);
6991 void theta_method_dispatcher::next_real_iter
6992 (
const model &md,
size_type ib,
const model::varnamelist &vl,
6993 const model::varnamelist &dl, model::real_matlist &matl,
6994 std::vector<model::real_veclist> &vectl,
6995 std::vector<model::real_veclist> &vectl_sym,
bool first_iter)
const {
6996 next_iter(md, ib, vl, dl, matl, vectl, vectl_sym, first_iter);
6999 void theta_method_dispatcher::next_complex_iter
7000 (
const model &md,
size_type ib,
const model::varnamelist &vl,
7001 const model::varnamelist &dl,
7002 model::complex_matlist &matl,
7003 std::vector<model::complex_veclist> &vectl,
7004 std::vector<model::complex_veclist> &vectl_sym,
7005 bool first_iter)
const {
7006 next_iter(md, ib, vl, dl, matl, vectl, vectl_sym, first_iter);
7009 void theta_method_dispatcher::asm_real_tangent_terms
7010 (
const model &md,
size_type ib, model::real_matlist &,
7011 std::vector<model::real_veclist> &,
7012 std::vector<model::real_veclist> &,
7013 build_version version)
const
7014 { md.brick_call(ib, version, 0); }
7016 void theta_method_dispatcher::asm_complex_tangent_terms
7017 (
const model &md,
size_type ib, model::complex_matlist &,
7018 std::vector<model::complex_veclist> &,
7019 std::vector<model::complex_veclist> &,
7020 build_version version)
const
7021 { md.brick_call(ib, version, 0); }
7023 theta_method_dispatcher::theta_method_dispatcher(
const std::string &THETA)
7024 : virtual_dispatcher(2) {
7025 param_names.push_back(THETA);
7029 (
model &md, dal::bit_vector ibricks,
const std::string &THETA) {
7030 pdispatcher pdispatch = std::make_shared<theta_method_dispatcher>(THETA);
7031 for (dal::bv_visitor i(ibricks); !i.finished(); ++i)
7036 (
model &md,
const std::string &U,
const std::string &V,
7037 const std::string &pdt,
const std::string &ptheta) {
7043 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for time step");
7045 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for parameter theta");
7048 scalar_type(1) - scalar_type(1) / theta[0]),
7051 scalar_type(1) / (theta[0]*dt[0])),
7054 -scalar_type(1) / (theta[0]*dt[0])),
7058 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for time step");
7059 const model_real_plain_vector &theta = md.
real_variable(ptheta);
7060 GMM_ASSERT1(gmm::vect_size(dt) == 1,
"Bad format for parameter theta");
7063 scalar_type(1) - scalar_type(1) / theta[0]),
7066 scalar_type(1) / (theta[0]*dt[0])),
7069 -scalar_type(1) / (theta[0]*dt[0])),
7082 const std::string &pdt,
const std::string &ptwobeta,
7083 const std::string &pgamma) {
7093 if (twobeta != gamma) {
7095 md.set_dispatch_coeff();
7099 md.
assembly(model::BUILD_RHS_WITH_LIN);
7102 model_complex_plain_vector W(nbdof), RHS(nbdof);
7103 gmm::copy(gmm::sub_vector(md.
complex_rhs(), md.interval_of_variable(U)),
7108 gmm::cg(md.linear_complex_matrix_term(id2dt2b, 0),
7109 W, RHS, gmm::identity_matrix(), iter);
7110 GMM_ASSERT1(iter.converged(),
"Velocity not well computed");
7112 gmm::scaled(W, complex_type(1)/(twobeta*dt)),
7116 if (twobeta != gamma) {
7118 md.set_dispatch_coeff();
7122 GMM_ASSERT1(
false,
"to be done");
7131 if (twobeta != gamma) {
7133 md.set_dispatch_coeff();
7137 md.
assembly(model::BUILD_RHS_WITH_LIN);
7140 model_real_plain_vector W(nbdof), RHS(nbdof);
7141 gmm::copy(gmm::sub_vector(md.
real_rhs(), md.interval_of_variable(U)),
7146 gmm::cg(md.linear_real_matrix_term(id2dt2b, 0),
7147 W, RHS, gmm::identity_matrix(), iter);
7148 GMM_ASSERT1(iter.converged(),
"Velocity not well computed");
7150 gmm::scaled(W, scalar_type(1)/(twobeta*dt)),
7154 if (twobeta != gamma) {
7156 md.set_dispatch_coeff();
7171 class midpoint_dispatcher :
public virtual_dispatcher {
7173 gmm::uint64_type id_num;
7177 typedef model::build_version build_version;
7179 void set_dispatch_coeff(
const model &md,
size_type ib)
const {
7180 md.matrix_coeff_of_brick(ib) = scalar_type(1)/scalar_type(2);
7181 md.rhs_coeffs_of_brick(ib)[0] = scalar_type(1);
7182 md.rhs_coeffs_of_brick(ib)[1] = scalar_type(1)/scalar_type(2);
7185 template <
typename MATLIST,
typename VECTLIST>
7186 inline void next_iter(
const model &md,
size_type ib,
7187 const model::varnamelist &vl,
7188 const model::varnamelist &dl,
7190 VECTLIST &vectl, VECTLIST &vectl_sym,
7191 bool first_iter)
const {
7193 pbrick pbr = md.brick_pointer(ib);
7197 if (!(pbr->is_linear()))
7198 md.add_temporaries(vl, id_num);
7199 md.add_temporaries(dl, id_num);
7204 if (pbr->is_linear()) {
7207 if (first_iter) md.update_brick(ib, model::BUILD_RHS);
7210 md.linear_brick_add_to_rhs(ib, 1, 0);
7215 (
const model &md,
size_type ib,
const model::varnamelist &vl,
7216 const model::varnamelist &dl, model::real_matlist &matl,
7217 std::vector<model::real_veclist> &vectl,
7218 std::vector<model::real_veclist> &vectl_sym,
bool first_iter)
const {
7219 next_iter(md, ib, vl, dl, matl, vectl, vectl_sym, first_iter);
7222 void next_complex_iter
7223 (
const model &md,
size_type ib,
const model::varnamelist &vl,
7224 const model::varnamelist &dl,
7225 model::complex_matlist &matl,
7226 std::vector<model::complex_veclist> &vectl,
7227 std::vector<model::complex_veclist> &vectl_sym,
7228 bool first_iter)
const {
7229 next_iter(md, ib, vl, dl, matl, vectl, vectl_sym, first_iter);
7232 void asm_real_tangent_terms
7233 (
const model &md,
size_type ib, model::real_matlist &,
7234 std::vector<model::real_veclist> &vectl,
7235 std::vector<model::real_veclist> &vectl_sym,
7236 build_version version)
const {
7238 scalar_type half = scalar_type(1)/scalar_type(2);
7239 pbrick pbr = md.brick_pointer(ib);
7242 const model::varnamelist &vl = md.varnamelist_of_brick(ib);
7243 const model::varnamelist &dl = md.datanamelist_of_brick(ib);
7245 if (!(pbr->is_linear())) {
7246 for (
size_type i = 0; i < vl.size(); ++i) {
7247 bool is_uptodate = md.temporary_uptodate(vl[i], id_num, ind);
7248 if (!is_uptodate && ind !=
size_type(-1))
7249 gmm::add(gmm::scaled(md.real_variable(vl[i], 0), half),
7250 gmm::scaled(md.real_variable(vl[i], 1), half),
7251 md.set_real_variable(vl[i], ind));
7252 md.set_default_iter_of_variable(vl[i], ind);
7257 for (
size_type i = 0; i < dl.size(); ++i) {
7258 bool is_uptodate = md.temporary_uptodate(dl[i], id_num, ind);
7259 if (!is_uptodate && ind !=
size_type(-1)) {
7260 gmm::add(gmm::scaled(md.real_variable(dl[i], 0), half),
7261 gmm::scaled(md.real_variable(dl[i], 1), half),
7262 md.set_real_variable(dl[i], ind));
7264 md.set_default_iter_of_variable(dl[i], ind);
7268 md.brick_call(ib, version, 0);
7269 if (pbr->is_linear()) {
7273 md.linear_brick_add_to_rhs(ib, 1, 1);
7276 md.reset_default_iter_of_variables(dl);
7277 if (!(pbr->is_linear()))
7278 md.reset_default_iter_of_variables(vl);
7281 virtual void asm_complex_tangent_terms
7282 (
const model &md,
size_type ib, model::complex_matlist &,
7283 std::vector<model::complex_veclist> &vectl,
7284 std::vector<model::complex_veclist> &vectl_sym,
7285 build_version version)
const {
7287 scalar_type half = scalar_type(1)/scalar_type(2);
7288 pbrick pbr = md.brick_pointer(ib);
7291 const model::varnamelist &vl = md.varnamelist_of_brick(ib);
7292 const model::varnamelist &dl = md.datanamelist_of_brick(ib);
7294 if (!(pbr->is_linear())) {
7295 for (
size_type i = 0; i < vl.size(); ++i) {
7296 bool is_uptodate = md.temporary_uptodate(vl[i], id_num, ind);
7297 if (!is_uptodate && ind !=
size_type(-1))
7298 gmm::add(gmm::scaled(md.complex_variable(vl[i], 0), half),
7299 gmm::scaled(md.complex_variable(vl[i], 1), half),
7300 md.set_complex_variable(vl[i], ind));
7301 md.set_default_iter_of_variable(vl[i], ind);
7306 for (
size_type i = 0; i < dl.size(); ++i) {
7307 bool is_uptodate = md.temporary_uptodate(dl[i], id_num, ind);
7308 if (!is_uptodate && ind !=
size_type(-1)) {
7309 gmm::add(gmm::scaled(md.complex_variable(dl[i], 0), half),
7310 gmm::scaled(md.complex_variable(dl[i], 1), half),
7311 md.set_complex_variable(dl[i], ind));
7313 md.set_default_iter_of_variable(dl[i], ind);
7317 md.brick_call(ib, version, 0);
7318 if (pbr->is_linear()) {
7322 md.linear_brick_add_to_rhs(ib, 1, 1);
7325 md.reset_default_iter_of_variables(dl);
7326 if (!(pbr->is_linear()))
7327 md.reset_default_iter_of_variables(vl);
7330 midpoint_dispatcher() : virtual_dispatcher(2)
7331 { id_num = act_counter(); }
7336 pdispatcher pdispatch = std::make_shared<midpoint_dispatcher>();
7337 for (dal::bv_visitor i(ibricks); !i.finished(); ++i)