GetFEM  5.4.2
getfem_models.cc
1 /*===========================================================================
2 
3  Copyright (C) 2009-2020 Yves Renard
4 
5  This file is a part of GetFEM
6 
7  GetFEM is free software; you can redistribute it and/or modify it
8  under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version along with the GCC Runtime Library
11  Exception either version 3.1 or (at your option) any later version.
12  This program is distributed in the hope that it will be useful, but
13  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  License and GCC Runtime Library Exception for more details.
16  You should have received a copy of the GNU Lesser General Public License
17  along with this program; if not, write to the Free Software Foundation,
18  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
19 
20 ===========================================================================*/
21 
22 #include <iomanip>
23 #include "gmm/gmm_range_basis.h"
24 #include "gmm/gmm_solver_cg.h"
26 #include "getfem/getfem_models.h"
33 
34 
35 namespace getfem {
36 
37  model::model(bool comp_version) {
38  init(); complex_version = comp_version;
39  is_linear_ = is_symmetric_ = is_coercive_ = true;
40  leading_dim = 0;
41  time_integration = 0; init_step = false; time_step = scalar_type(1);
45  ("neighbor_element", interpolate_transformation_neighbor_instance());
46 
47  ga_tree tree1;
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;
53  tree1.root->nbc1 = 0;
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);
58 
59  ga_tree tree2;
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;
65  tree2.root->nbc1 = 0;
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);
70  }
71 
72  void model::var_description::set_size() {
73  clear_temporaries();
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()
79  : 1);
80  s *= qdim();
81  for (size_type i = 0; i < n_iter; ++i)
82  if (is_complex)
83  complex_value[i].resize(s);
84  else
85  real_value[i].resize(s);
86  if (is_affine_dependent) {
87  if (is_complex)
88  affine_complex_value.resize(s);
89  else
90  affine_real_value.resize(s);
91  }
92  }
93 
94  size_type model::var_description::add_temporary(gmm::uint64_type id_num) {
95  size_type nit = n_iter;
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) {
99  ++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);
103  v_num_iter[nit] = 0;
104  if (is_complex) {
105  size_type s = complex_value[0].size();
106  complex_value.resize(n_iter + n_temp_iter);
107  complex_value[nit].resize(s);
108  } else {
109  size_type s = real_value[0].size();
110  real_value.resize(n_iter + n_temp_iter);
111  real_value[nit].resize(s);
112  }
113  }
114  return nit;
115  }
116 
117  void model::var_description::clear_temporaries() {
118  n_temp_iter = 0;
119  default_iter = 0;
120  if (is_complex)
121  complex_value.resize(n_iter);
122  else
123  real_value.resize(n_iter);
124  }
125 
126  bool model::check_name_validity(const std::string &name, bool assert) const {
127 
128  if (variables.count(name) != 0) {
129  GMM_ASSERT1(!assert, "Variable " << name << " already exists");
130  return false;
131  } else if (variable_groups.count(name) != 0) {
132  GMM_ASSERT1(!assert,
133  name << " corresponds to an already existing group of "
134  "variables name");
135  return false;
136  } else if (macro_exists(name)) {
137  GMM_ASSERT1(!assert,
138  name << " corresponds to an already existing macro");
139  return false;
140  } else if (name.compare("X") == 0) {
141  GMM_ASSERT1(!assert, "X is a reserved keyword of the generic "
142  "assembly language");
143  return false;
144  }
145 
146  int ga_valid = ga_check_name_validity(name);
147  if (ga_valid == 1) {
148  GMM_ASSERT1(!assert, "Invalid variable name, corresponds to an "
149  "operator or function name of the generic assembly language");
150  return false;
151  } else if (ga_valid == 2) {
152  GMM_ASSERT1(!assert, "Invalid variable name having a reserved "
153  "prefix used by the generic assembly language");
154  return false;
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()) {
159  GMM_ASSERT1(!assert,
160  "Dot and Previous are reserved prefix used for time "
161  "integration schemes");
162  return false;
163  }
164  }
165 
166  bool valid = !name.empty() && isalpha(name[0]);
167  if (valid)
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 << "\"");
172  return valid;
173  }
174 
175  std::string model::new_name(const std::string &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) {
179  std::stringstream m;
180  m << name << '_' << i;
181  res_name = m.str();
182  valid = check_name_validity(res_name, false);
183  }
184  for (size_type i = 2; !valid && i < 1000; ++i) {
185  std::stringstream m;
186  m << "new_" << name << '_' << i;
187  res_name = m.str();
188  valid = check_name_validity(res_name, false);
189  }
190  GMM_ASSERT1(valid, "Illegal variable name: " << name);
191  return res_name;
192  }
193 
194 
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);
199  return it;
200  }
201 
202  const model::var_description &
203  model::variable_description(const std::string &name) const {
204  return find_variable(name)->second;
205  }
206 
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);
210  }
211  if (!(v.compare(0, 3, "Dot")) && (v[3] == '_' || v[4] == '_')) {
212  v = v.substr((v[3] == '_') ? 4 : 5);
213  }
214  if (is_old(v)) v = no_old_prefix_name(v);
215  return v;
216  }
217 
218  bool is_old(const std::string &name){
219  return name.substr(0, PREFIX_OLD_LENGTH) == PREFIX_OLD;
220  }
221 
222  std::string no_old_prefix_name(const std::string &name){
223  return is_old(name) ? name.substr(PREFIX_OLD_LENGTH) : name;
224  }
225 
226  bool model::is_disabled_variable(const std::string &name) const {
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;
233  }
234 
235  bool model::is_data(const std::string &name) const {
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;
241  }
242 
243  bool model::is_true_data(const std::string &name) const {
244  return is_old(name) || !(variable_description(name).is_variable);
245  }
246 
247  bool model::is_internal_variable(const std::string &name) const {
248  if (is_old(name)) return false;
249  const auto &var_descr = variable_description(name);
250  return var_descr.is_internal && var_descr.is_enabled();
251  }
252 
253  bool model::is_affine_dependent_variable(const std::string &name) const {
254  return !(is_old(name)) && variable_description(name).is_affine_dependent;
255  }
256 
257  const std::string &
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;
262  }
263 
264  const scalar_type &
265  model::factor_of_variable(const std::string &name) const {
266  return variable_description(name).alpha;
267  }
268 
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();
275  }
276  }
277 
278  bool model::is_im_data(const std::string &name) const {
279  return variable_description(no_old_prefix_name(name)).imd != 0;
280  }
281 
282  const im_data *
283  model::pim_data_of_variable(const std::string &name) const {
284  return variable_description(no_old_prefix_name(name)).imd;
285  }
286 
287  const gmm::uint64_type &
288  model::version_number_of_data_variable(const std::string &name,
289  size_type niter) const {
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];
293  }
294 
295  size_type model::nb_dof(bool with_internal) const {
296  context_check();
297  if (act_size_to_be_done) actualize_sizes();
298  if (complex_version)
299  return gmm::vect_size(crhs);
300  else if (with_internal && gmm::vect_size(full_rrhs))
301  return gmm::vect_size(full_rrhs);
302  else
303  return gmm::vect_size(rrhs);
304  }
305 
306  void model::resize_global_system() const {
307 
308  size_type full_size = 0;
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();
316  }
317  }
318  size_type primary_size = full_size;
319 
320  for (auto &&v : variables)
321  if (v.second.is_internal && v.second.is_enabled()) { // is_internal_variable()
322  v.second.I = gmm::sub_interval(full_size, v.second.size());
323  full_size += v.second.size();
324  }
325 
326  for (auto &&v : variables)
327  if (v.second.is_affine_dependent) {
328  v.second.I = variables.find(v.second.org_name)->second.I;
329  v.second.set_size();
330  }
331 
332  if (complex_version) {
333  gmm::resize(cTM, primary_size, primary_size);
334  gmm::resize(crhs, primary_size);
335  }
336  else {
337  gmm::resize(rTM, primary_size, primary_size);
338  gmm::resize(rrhs, primary_size);
339  }
340 
341  if (full_size > primary_size) {
342  GMM_ASSERT1(has_internal_variables(), "Internal error");
343  gmm::resize(internal_rTM, full_size-primary_size, primary_size);
344  gmm::resize(full_rrhs, full_size);
345  gmm::resize(internal_sol, full_size-primary_size);
346  } else {
347  GMM_ASSERT1(!(has_internal_variables()), "Internal error");
348  gmm::resize(internal_rTM, 0, 0);
349  full_rrhs.clear();
350  }
351 
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;
356  break;
357  }
358  }
359 
360  void model::actualize_sizes() const {
361  // cout << "begin act size" << endl;
362  bool actualized = false;
363  getfem::local_guard lock = locks_.get_lock();
364  if (actualized) return; // If multiple threads are calling the method
365 
366  act_size_to_be_done = false;
367 
368  std::map<std::string, std::vector<std::string> > multipliers;
369  std::set<std::string> tobedone;
370 
371 // #if GETFEM_PARA_LEVEL > 1
372 // int rk; MPI_Comm_rank(MPI_COMM_WORLD, &rk);
373 // double t_ref = MPI_Wtime();
374 // cout << "Actualize size called from thread " << rk << endl;
375 // #endif
376 
377 
378  // In case of change in fems or mims, linear terms have to be recomputed
379  // We could select which brick is to be recomputed if we would be able
380  // to know which fem or mim is changed
381  // -> already taken into account in update_brick()
382  // for (dal::bv_visitor ib(valid_bricks); !ib.finished(); ++ib)
383  // bricks[ib].terms_to_be_computed = true;
384 
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);
400  }
401  }
402  switch (vdescr.filter) {
403  case VDESCRFILTER_NO:
404  if (vdescr.v_num < vdescr.mf->version_number()) {
405  vdescr.set_size();
406  vdescr.v_num = act_counter();
407  }
408  break;
409  case VDESCRFILTER_REGION:
410  if (vdescr.v_num < vdescr.mf->version_number()) {
411  dal::bit_vector
412  dor = vdescr.mf->dof_on_region(vdescr.filter_region);
413  vdescr.partial_mf->adapt(dor);
414  vdescr.set_size();
415  vdescr.v_num = act_counter();
416  }
417  break;
418  default : break;
419  }
420  }
421 
422  if (vdescr.imd != 0
423  && vdescr.v_num < vdescr.imd->version_number()) {
424  vdescr.set_size();
425  vdescr.v_num = act_counter();
426  }
427  }
428 
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)) {
435  // This step forces the recomputation of corresponding bricks.
436  // A test to check if a modification is really necessary could
437  // be done first ... (difficult to coordinate with other
438  // multipliers)
439  dal::bit_vector alldof; alldof.add(0, vdescr.mf->nb_dof());
440  vdescr.partial_mf->adapt(alldof);
441  vdescr.set_size();
442  vdescr.v_num = act_counter();
443  }
444  }
445  }
446 
447  resize_global_system();
448 
449  for (const std::string &vname : tobedone) {
450 // #if GETFEM_PARA_LEVEL > 1
451 // double tt_ref = MPI_Wtime();
452 // if (!rk) cout << "compute size of multipliers for " << vname
453 // << endl;
454 // #endif
455 
456  const std::vector<std::string> &mults = multipliers[vname];
457  const var_description &vdescr = variable_description(vname);
458 
459  gmm::col_matrix< gmm::rsvector<scalar_type> > MGLOB;
460  if (mults.size() > 1) {
461  size_type s = 0;
462  for (const std::string &mult : mults)
463  s += variable_description(mult).mf->nb_dof();
464  gmm::resize(MGLOB, vdescr.mf->nb_dof(), s);
465  }
466  size_type s = 0;
467  std::set<size_type> glob_columns;
468  for (const std::string &multname : mults) {
469  var_description &multdescr = variables.find(multname)->second;
470 
471  // Obtaining the coupling matrix between the multipier and
472  // the primal variable. A search is done on all the terms of the
473  // model. Only the corresponding linear terms are added.
474  // If no linear term is available, a mass matrix is used.
475  gmm::col_matrix< gmm::rsvector<scalar_type> >
476  MM(vdescr.associated_mf().nb_dof(), multdescr.mf->nb_dof());
477  bool termadded = false;
478 
479  if (multdescr.filter & VDESCRFILTER_CTERM) {
480 
481  for (dal::bv_visitor ib(valid_bricks); !ib.finished(); ++ib) {
482  const brick_description &brick = bricks[ib];
483  bool bupd = false;
484  bool cplx = is_complex() && brick.pbr->is_complex();
485 
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;
491  }
492  if (multc && varc) {
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");
499  {
500  gmm::clear(rTM);
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);
509  );
510  } //accumulated_distro scope
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);
514  bupd = false;
515  }
516  }
517  }
518 
519  for (size_type j = 0; j < brick.tlist.size(); ++j) {
520  const term_description &term = brick.tlist[j];
521 
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;
528  }
529  if (multc && varc) {
530  GMM_ASSERT1(!cplx, "Sorry, not taken into account");
531  generic_expressions.clear();
532  if (!bupd) {
533  brick.terms_to_be_computed = true;
534  update_brick(ib, BUILD_MATRIX);
535  bupd = true;
536  }
537  gmm::add(gmm::sub_matrix(brick.rmatlist[j],
538  multdescr.I, vdescr.I),
539  MM);
540  gmm::add(gmm::transposed(gmm::sub_matrix
541  (brick.rmatlist[j],
542  vdescr.I, multdescr.I)),
543  MM);
544  termadded = true;
545  }
546  } else if (multname.compare(term.var1) == 0 &&
547  vname.compare(term.var2) == 0) {
548  if (!bupd) {
549  brick.terms_to_be_computed = true;
550  update_brick(ib, BUILD_MATRIX);
551  bupd = true;
552  }
553  if (cplx)
554  gmm::add(gmm::transposed(gmm::real_part(brick.cmatlist[j])),
555  MM);
556  else
557  gmm::add(gmm::transposed(brick.rmatlist[j]), MM);
558  termadded = true;
559 
560  } else if (multname.compare(term.var2) == 0 &&
561  vname.compare(term.var1) == 0) {
562  if (!bupd) {
563  brick.terms_to_be_computed = true;
564  update_brick(ib, BUILD_MATRIX);
565  bupd = true;
566  }
567  if (cplx)
568  gmm::add(gmm::real_part(brick.cmatlist[j]), MM);
569  else
570  gmm::add(brick.rmatlist[j], MM);
571  termadded = true;
572  }
573  }
574  }
575  }
576 
577  if (!termadded)
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);
583  asm_mass_matrix(MM, *(multdescr.filter_mim), vdescr.associated_mf(),
584  *(multdescr.mf), rg);
585  }
586 
587  MPI_SUM_SPARSE_MATRIX(MM);
588 
589  //
590  // filtering
591  //
592  std::set<size_type> columns;
593  gmm::range_basis(MM, columns);
594  if (columns.size() == 0)
595  GMM_WARNING1("Empty basis found for multiplier " << multname);
596 
597  if (mults.size() > 1) {
598  gmm::copy(MM, gmm::sub_matrix
599  (MGLOB,
600  gmm::sub_interval(0, vdescr.associated_mf().nb_dof()),
601  gmm::sub_interval(s, multdescr.mf->nb_dof())));
602  for (const size_type &icol : columns)
603  glob_columns.insert(s + icol);
604  s += multdescr.mf->nb_dof();
605  } else {
606  dal::bit_vector kept;
607  for (const size_type &icol : columns)
608  kept.add(icol);
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();
614  }
615  }
616 
617 // #if GETFEM_PARA_LEVEL > 1
618 // if (!rk) cout << "Range basis for multipliers for " << vname << " time : " << MPI_Wtime()-tt_ref << endl;
619 // #endif
620 
621  if (mults.size() > 1) {
622  gmm::range_basis(MGLOB, glob_columns, 1E-12, gmm::col_major(), true);
623 
624 // #if GETFEM_PARA_LEVEL > 1
625 // if (!rk) cout << "Producing partial mf for multipliers for " << vname << " time : " << MPI_Wtime()-tt_ref << endl;
626 // #endif
627 
628  s = 0;
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();
641  }
642  }
643 // #if GETFEM_PARA_LEVEL > 1
644 // if (!rk) cout << "End compute size of multipliers for " << vname << " time : " << MPI_Wtime()-tt_ref << endl;
645 // #endif
646  }
647 
648  resize_global_system();
649  actualized = true;
650 // #if GETFEM_PARA_LEVEL > 1
651 // cout << "Actualize sizes time from thread " << rk << " : " << MPI_Wtime()-t_ref << endl;
652 
653 // #endif
654 
655  // cout << "end act size" << endl;
656  }
657 
658 
659  void model::listvar(std::ostream &ost) const {
660  if (variables.size() == 0)
661  ost << "Model with no variable nor data" << endl;
662  else {
663  ost << "List of model variables and data:" << endl;
664  for (int vartype=0; vartype < 3; ++vartype)
665  for (const auto &v : variables) {
666  const var_description &vdescr = v.second;
667  bool is_variable = vdescr.is_variable;
668  bool is_disabled = is_variable && is_disabled_variable(v.first);
669  if (vartype == 0) { // Only enabled variables
670  if (!is_variable || is_disabled) continue;
671  } else if (vartype == 1) { // Only disabled variables
672  if (!is_disabled) continue;
673  } else if (vartype == 2) { // Only data
674  if (is_variable) continue;
675  }
676  ost << (is_variable ? "Variable " : "Data ");
677  ost << std::setw(30) << std::left << v.first;
678  ost << std::setw(2) << std::right << vdescr.n_iter;
679  ost << ((vdescr.n_iter == 1) ? " copy " : " copies ");
680  ost << (vdescr.is_fem_dofs ? "fem dependant " : "constant size ");
681  ost << std::setw(8) << std::right << vdescr.size();
682  if (is_complex()) ost << " complex";
683  ost << ((vdescr.size() > 1) ? " doubles." : " double.");
684  ost << (is_disabled ? "\t (disabled)" : "\t ");
685  if (vdescr.imd != 0) ost << "\t (is im_data)";
686  if (vdescr.is_affine_dependent) ost << "\t (is affine dependent)";
687  ost << endl;
688  }
689  for (const auto &vargroup : variable_groups) {
690  ost << "Variable group " << std::setw(30) << std::left
691  << vargroup.first;
692  if (vargroup.second.size()) {
693  bool first(true);
694  for (const std::string &vname : vargroup.second) {
695  ost << (first ? " " : ", ") << vname;
696  first = false;
697  }
698  ost << endl;
699  } else
700  ost << " empty" << endl;
701  }
702  }
703  }
704 
705  void model::listresiduals(std::ostream &ost) const {
706  context_check(); if (act_size_to_be_done) actualize_sizes();
707  if (variables.size() == 0)
708  ost << "Model with no variable nor data" << endl;
709  else {
710  bool firstvar(true);
711  for (const auto &v : variables) {
712  if (v.second.is_variable) {
713  const model_real_plain_vector &rhs = v.second.is_internal
714  ? full_rrhs : rrhs;
715  const gmm::sub_interval &II = interval_of_variable(v.first);
716  scalar_type res = gmm::vect_norm2(gmm::sub_vector(rhs, II));
717  if (!firstvar) cout << ", ";
718  ost << "res_" << v.first << "= " << std::setw(11) << res;
719  firstvar = false;
720  }
721  }
722  ost << endl;
723  }
724  }
725 
726  void model::add_fixed_size_variable(const std::string &name, size_type size,
727  size_type niter) {
728  bgeot::multi_index sizes(1);
729  sizes[0] = size;
730  add_fixed_size_variable(name, sizes, niter);
731  }
732 
733  void model::add_fixed_size_variable(const std::string &name,
734  const bgeot::multi_index &sizes,
735  size_type niter) {
736  check_name_validity(name);
737  variables.emplace(name, var_description(true, is_complex(), 0, 0, niter));
738  variables[name].qdims = sizes;
739  act_size_to_be_done = true;
740  variables[name].set_size();
741  GMM_ASSERT1(variables[name].qdim(),
742  "Variables of null size are not allowed");
743  }
744 
745  void model::resize_fixed_size_variable(const std::string &name,
746  size_type size) {
747  bgeot::multi_index sizes(1);
748  sizes[0] = size;
749  resize_fixed_size_variable(name, sizes);
750  }
751 
752  void model::resize_fixed_size_variable(const std::string &name,
753  const bgeot::multi_index &sizes) {
754  GMM_ASSERT1(!(variables[name].is_fem_dofs),
755  "Cannot explicitly resize a fem variable or data");
756  GMM_ASSERT1(variables[name].imd == 0,
757  "Cannot explicitly resize an im data");
758  variables[name].qdims = sizes;
759  variables[name].set_size();
760  }
761 
762  void model::add_fixed_size_data(const std::string &name, size_type size,
763  size_type niter) {
764  bgeot::multi_index sizes(1);
765  sizes[0] = size;
766  add_fixed_size_data(name, sizes, niter);
767  }
768 
769  void model::add_fixed_size_data(const std::string &name,
770  const bgeot::multi_index &sizes,
771  size_type niter) {
772  check_name_validity(name);
773  variables.emplace(name, var_description(false, is_complex(), 0, 0, niter));
774  variables[name].qdims = sizes;
775  GMM_ASSERT1(variables[name].qdim(), "Data of null size are not allowed");
776  variables[name].set_size();
777  }
778 
779  void model::add_initialized_matrix_data(const std::string &name,
780  const base_matrix &M) {
781  add_fixed_size_data(name, bgeot::multi_index(gmm::mat_nrows(M),
782  gmm::mat_ncols(M)));
783  GMM_ASSERT1(!(is_complex()), "Sorry, complex version to be done");
784  gmm::copy(M.as_vector(), set_real_variable(name));
785  }
786 
787  void model::add_initialized_matrix_data(const std::string &name,
788  const base_complex_matrix &M) {
789  add_fixed_size_data(name, bgeot::multi_index(gmm::mat_nrows(M),
790  gmm::mat_ncols(M)));
791  GMM_ASSERT1(!(is_complex()), "Sorry, complex version to be done");
792  gmm::copy(M.as_vector(), set_complex_variable(name));
793  }
794 
795  void model::add_initialized_tensor_data(const std::string &name,
796  const base_tensor &t) {
797  add_fixed_size_data(name, t.sizes(), 1);
798  GMM_ASSERT1(!(is_complex()), "Sorry, complex version to be done");
799  gmm::copy(t.as_vector(), set_real_variable(name));
800  }
801 
802  void model::add_initialized_tensor_data(const std::string &name,
803  const base_complex_tensor &t) {
804  add_fixed_size_data(name, t.sizes(), 1);
805  GMM_ASSERT1(!(is_complex()), "Sorry, complex version to be done");
806  gmm::copy(t.as_vector(), set_complex_variable(name));
807  }
808 
809  void model::add_im_variable(const std::string &name, const im_data &imd,
810  size_type niter) {
811  check_name_validity(name);
812  variables.emplace(name,
813  var_description(true, is_complex(), 0, &imd, niter));
814  variables[name].set_size();
815  add_dependency(imd);
816  act_size_to_be_done = true;
817  }
818 
819  void model::add_internal_im_variable(const std::string &name,
820  const im_data &imd) {
821  add_im_variable(name, imd);
822  variables[name].is_internal = true;
823  }
824 
825  void model::add_im_data(const std::string &name, const im_data &imd,
826  size_type niter) {
827  check_name_validity(name);
828  variables.emplace(name,
829  var_description(false, is_complex(), 0, &imd, niter));
830  variables[name].set_size();
831  add_dependency(imd);
832  }
833 
834  void model::add_fem_variable(const std::string &name, const mesh_fem &mf,
835  size_type niter) {
836  check_name_validity(name);
837  variables.emplace(name, var_description(true, is_complex(), &mf, 0, niter,
838  VDESCRFILTER_NO));
839  variables[name].set_size();
840  add_dependency(mf);
841  act_size_to_be_done = true;
842  leading_dim = std::max(leading_dim, mf.linked_mesh().dim());
843  }
844 
845  void model::add_filtered_fem_variable(const std::string &name,
846  const mesh_fem &mf,
847  size_type region, size_type niter) {
848  check_name_validity(name);
849  variables.emplace(name, var_description(true, is_complex(), &mf, 0, niter,
850  VDESCRFILTER_REGION, region));
851  variables[name].set_size();
852  act_size_to_be_done = true;
853  add_dependency(mf);
854  }
855 
856  void model::add_affine_dependent_variable(const std::string &name,
857  const std::string &org_name,
858  scalar_type alpha) {
859  check_name_validity(name);
860  VAR_SET::const_iterator it = find_variable(org_name);
861  GMM_ASSERT1(it->second.is_variable && !(it->second.is_affine_dependent),
862  "The original variable should be a variable");
863  variables.emplace(name, variables[org_name]);
864  variables[name].is_affine_dependent = true;
865  variables[name].org_name = org_name;
866  variables[name].alpha = alpha;
867  variables[name].set_size();
868  }
869 
870  void model::add_fem_data(const std::string &name, const mesh_fem &mf,
871  dim_type qdim, size_type niter) {
872  bgeot::multi_index sizes(1);
873  sizes[0] = qdim;
874  add_fem_data(name, mf, sizes, niter);
875  }
876 
877  void model::add_fem_data(const std::string &name, const mesh_fem &mf,
878  const bgeot::multi_index &sizes, size_type niter) {
879  check_name_validity(name);
880  variables.emplace(name, var_description(false, is_complex(), &mf, 0, niter,
881  VDESCRFILTER_NO));
882  variables[name].qdims = sizes;
883  GMM_ASSERT1(variables[name].qdim(), "Data of null size are not allowed");
884  variables[name].set_size();
885  add_dependency(mf);
886  }
887 
888  void model::add_multiplier(const std::string &name, const mesh_fem &mf,
889  const std::string &primal_name,
890  size_type niter) {
891  check_name_validity(name);
892  variables.emplace(name, var_description(true, is_complex(), &mf, 0, niter,
893  VDESCRFILTER_CTERM, size_type(-1),
894  primal_name));
895  variables[name].set_size();
896  act_size_to_be_done = true;
897  add_dependency(mf);
898  }
899 
900  void model::add_multiplier(const std::string &name, const mesh_fem &mf,
901  size_type region, const std::string &primal_name,
902  size_type niter) {
903  check_name_validity(name);
904  variables.emplace(name, var_description(true, is_complex(), &mf, 0, niter,
905  VDESCRFILTER_REGION_CTERM, region,
906  primal_name));
907  variables[name].set_size();
908  act_size_to_be_done = true;
909  add_dependency(mf);
910  }
911 
912  void model::add_multiplier(const std::string &name, const mesh_fem &mf,
913  const std::string &primal_name,
914  const mesh_im &mim,
915  size_type region, size_type niter) {
916  check_name_validity(name);
917  variables.emplace(name, var_description(true, is_complex(), &mf, 0, niter,
918  VDESCRFILTER_INFSUP, region,
919  primal_name, &mim));
920  variables[name].set_size();
921  act_size_to_be_done = true;
922  add_dependency(mf);
923  add_dependency(mim);
924  }
925 
926  void model::disable_variable(const std::string &name) {
927  enable_variable(name, false);
928  }
929 
930  void model::enable_variable(const std::string &name, bool enabled) {
931  VAR_SET::iterator it = variables.find(name);
932  GMM_ASSERT1(it != variables.end(), "Undefined variable " << name);
933  it->second.is_disabled = !enabled;
934  for (auto &&v : variables) {
935  if (((v.second.filter & VDESCRFILTER_INFSUP) ||
936  (v.second.filter & VDESCRFILTER_CTERM))
937  && name.compare(v.second.filter_var) == 0) {
938  v.second.is_disabled = !enabled;
939  }
940  if (v.second.is_variable && v.second.is_affine_dependent
941  && name.compare(v.second.org_name) == 0)
942  v.second.is_disabled = !enabled;
943  }
944  if (!act_size_to_be_done) resize_global_system();
945  }
946 
947  bool model::variable_exists(const std::string &name) const {
948  return variables.count(no_old_prefix_name(name)) > 0;
949  }
950 
951  void model::add_macro(const std::string &name, const std::string &expr) {
952  check_name_validity(name.substr(0, name.find("(")));
953  macro_dict.add_macro(name, expr);
954  }
955 
956  void model::del_macro(const std::string &name)
957  { macro_dict.del_macro(name); }
958 
960  GMM_ASSERT1(valid_bricks[ib], "Inexistent brick");
961  valid_bricks.del(ib);
962  active_bricks.del(ib);
963 
964  for (size_type i = 0; i < bricks[ib].mims.size(); ++i) {
965  const mesh_im *mim = bricks[ib].mims[i];
966  bool found = false;
967  for (dal::bv_visitor ibb(valid_bricks); !ibb.finished(); ++ibb) {
968  for (size_type j = 0; j < bricks[ibb].mims.size(); ++j)
969  if (bricks[ibb].mims[j] == mim) found = true;
970  }
971  for (const auto &v : variables) {
972  if (v.second.is_fem_dofs &&
973  (v.second.filter & VDESCRFILTER_INFSUP) &&
974  mim == v.second.filter_mim) found = true;
975  }
976  if (!found) sup_dependency(*mim);
977  }
978 
979  is_linear_ = is_symmetric_ = is_coercive_ = true;
980  for (dal::bv_visitor ibb(valid_bricks); !ibb.finished(); ++ibb) {
981  is_linear_ = is_linear_ && bricks[ibb].pbr->is_linear();
982  is_symmetric_ = is_symmetric_ && bricks[ibb].pbr->is_symmetric();
983  is_coercive_ = is_coercive_ && bricks[ibb].pbr->is_coercive();
984  }
985  bricks[ib] = brick_description();
986  }
987 
988  void model::delete_variable(const std::string &varname) {
989  for (dal::bv_visitor ibb(valid_bricks); !ibb.finished(); ++ibb) {
990  for (const auto &vname : bricks[ibb].vlist)
991  GMM_ASSERT1(varname.compare(vname),
992  "Cannot delete a variable which is still used by a brick");
993  for (const auto &dname : bricks[ibb].dlist)
994  GMM_ASSERT1(varname.compare(dname),
995  "Cannot delete a data which is still used by a brick");
996  }
997 
998  VAR_SET::const_iterator it = find_variable(varname);
999 
1000  if (it->second.is_fem_dofs) {
1001  const mesh_fem *mf = it->second.mf;
1002  bool found = false;
1003  for(VAR_SET::iterator it2 = variables.begin();
1004  it2 != variables.end(); ++it2) {
1005  if (it != it2 && it2->second.is_fem_dofs && mf == it2->second.mf)
1006  found = true;
1007  }
1008  if (!found) sup_dependency(*mf);
1009 
1010  if (it->second.filter & VDESCRFILTER_INFSUP) {
1011  const mesh_im *mim = it->second.filter_mim;
1012  found = false;
1013  for (dal::bv_visitor ibb(valid_bricks); !ibb.finished(); ++ibb) {
1014  for (size_type j = 0; j < bricks[ibb].mims.size(); ++j)
1015  if (bricks[ibb].mims[j] == mim) found = true;
1016  }
1017  for (VAR_SET::iterator it2 = variables.begin();
1018  it2 != variables.end(); ++it2) {
1019  if (it != it2 && it2->second.is_fem_dofs &&
1020  (it2->second.filter & VDESCRFILTER_INFSUP) &&
1021  mim == it2->second.filter_mim) found = true;
1022  }
1023  if (!found) sup_dependency(*mim);
1024  }
1025  }
1026 
1027  if (it->second.imd != 0) sup_dependency(*(it->second.imd));
1028 
1029  variables.erase(varname);
1030  act_size_to_be_done = true;
1031  }
1032 
1033  size_type model::add_brick(pbrick pbr, const varnamelist &varnames,
1034  const varnamelist &datanames,
1035  const termlist &terms,
1036  const mimlist &mims, size_type region) {
1037  size_type ib = valid_bricks.first_false();
1038 
1039  for (size_type i = 0; i < terms.size(); ++i)
1040  if (terms[i].is_global && terms[i].is_matrix_term && pbr->is_linear())
1041  GMM_ASSERT1(false, "Global linear matrix terms are not allowed");
1042 
1043  if (ib == bricks.size())
1044  bricks.push_back(brick_description(pbr, varnames, datanames, terms,
1045  mims, region));
1046  else
1047  bricks[ib] = brick_description(pbr, varnames, datanames, terms,
1048  mims, region);
1049  active_bricks.add(ib);
1050  valid_bricks.add(ib);
1051 
1052  // The brick itself already reacts to a mesh_im change in update_brick()
1053  // for (size_type i = 0; i < bricks[ib].mims.size(); ++i)
1054  // add_dependency(*(bricks[ib].mims[i]));
1055 
1056  GMM_ASSERT1(pbr->is_real() || is_complex(),
1057  "Impossible to add a complex brick to a real model");
1058  if (is_complex() && pbr->is_complex()) {
1059  bricks[ib].cmatlist.resize(terms.size());
1060  bricks[ib].cveclist[0].resize(terms.size());
1061  bricks[ib].cveclist_sym[0].resize(terms.size());
1062  } else {
1063  bricks[ib].rmatlist.resize(terms.size());
1064  bricks[ib].rveclist[0].resize(terms.size());
1065  bricks[ib].rveclist_sym[0].resize(terms.size());
1066  }
1067  is_linear_ = is_linear_ && pbr->is_linear();
1068  is_symmetric_ = is_symmetric_ && pbr->is_symmetric();
1069  is_coercive_ = is_coercive_ && pbr->is_coercive();
1070 
1071  for (const auto &vname : varnames)
1072  GMM_ASSERT1(variables.count(vname),
1073  "Undefined model variable " << vname);
1074  // cout << "dl == " << datanames << endl;
1075  for (const auto &dname : datanames)
1076  GMM_ASSERT1(variables.count(dname),
1077  "Undefined model data or variable " << dname);
1078 
1079  return ib;
1080  }
1081 
1083  GMM_ASSERT1(valid_bricks[ib], "Inexistent brick");
1084  touch_brick(ib);
1085  bricks[ib].mims.push_back(&mim);
1086  add_dependency(mim);
1087  }
1088 
1089  void model::change_terms_of_brick(size_type ib, const termlist &terms) {
1090  GMM_ASSERT1(valid_bricks[ib], "Inexistent brick");
1091  touch_brick(ib);
1092  bricks[ib].tlist = terms;
1093  if (is_complex() && bricks[ib].pbr->is_complex()) {
1094  bricks.back().cmatlist.resize(terms.size());
1095  bricks.back().cveclist[0].resize(terms.size());
1096  bricks.back().cveclist_sym[0].resize(terms.size());
1097  } else {
1098  bricks.back().rmatlist.resize(terms.size());
1099  bricks.back().rveclist[0].resize(terms.size());
1100  bricks.back().rveclist_sym[0].resize(terms.size());
1101  }
1102  }
1103 
1104  void model::change_variables_of_brick(size_type ib, const varnamelist &vl) {
1105  GMM_ASSERT1(valid_bricks[ib], "Inexistent brick");
1106  touch_brick(ib);
1107  bricks[ib].vlist = vl;
1108  for (const auto &v : vl)
1109  GMM_ASSERT1(variables.count(v), "Undefined model variable " << v);
1110  }
1111 
1112  void model::change_data_of_brick(size_type ib, const varnamelist &dl) {
1113  GMM_ASSERT1(valid_bricks[ib], "Inexistent brick");
1114  touch_brick(ib);
1115  bricks[ib].dlist = dl;
1116  for (const auto &v : dl)
1117  GMM_ASSERT1(variables.count(v), "Undefined model variable " << v);
1118  }
1119 
1120  void model::change_mims_of_brick(size_type ib, const mimlist &ml) {
1121  GMM_ASSERT1(valid_bricks[ib], "Inexistent brick");
1122  touch_brick(ib);
1123  bricks[ib].mims = ml;
1124  for (const auto &mim : ml) add_dependency(*mim);
1125  }
1126 
1128  GMM_ASSERT1(valid_bricks[ib], "Inexistent brick");
1129  touch_brick(ib);
1130  bricks[ib].is_update_brick = flag;
1131  }
1132 
1133  void model::set_time(scalar_type t, bool to_init) {
1134  static const std::string varname("t");
1135  VAR_SET::iterator it = variables.find(varname);
1136  if (it == variables.end()) {
1137  add_fixed_size_data(varname, 1);
1138  } else {
1139  GMM_ASSERT1(it->second.size() == 1, "Time data should be of size 1");
1140  }
1141  if (it == variables.end() || to_init) {
1142  if (is_complex())
1143  set_complex_variable(varname)[0] = complex_type(t);
1144  else
1145  set_real_variable(varname)[0] = t;
1146  }
1147  }
1148 
1149  scalar_type model::get_time() {
1150  static const std::string varname("t");
1151  set_time(scalar_type(0), false);
1152  if (is_complex())
1153  return gmm::real(complex_variable(varname)[0]);
1154  else
1155  return real_variable(varname)[0];
1156  }
1157 
1158  void model::call_init_affine_dependent_variables(int version) {
1159  for (VAR_SET::iterator it = variables.begin();
1160  it != variables.end(); ++it) {
1161  var_description &vdescr = it->second;
1162  if (vdescr.is_variable && vdescr.ptsc) {
1163  if (version == 2)
1164  vdescr.ptsc->init_affine_dependent_variables_precomputation(*this);
1165  else
1166  vdescr.ptsc->init_affine_dependent_variables(*this);
1167  }
1168  }
1169  }
1170 
1171  void model::shift_variables_for_time_integration() {
1172  for (VAR_SET::iterator it = variables.begin();
1173  it != variables.end(); ++it)
1174  if (it->second.is_variable && it->second.ptsc)
1175  it->second.ptsc->shift_variables(*this);
1176  }
1177 
1178  void model::add_time_integration_scheme(const std::string &varname,
1179  ptime_scheme ptsc) {
1180  VAR_SET::iterator it = variables.find(varname);
1181  GMM_ASSERT1(it != variables.end(), "Undefined variable " << varname);
1182  GMM_ASSERT1(it->second.is_variable && !(it->second.is_affine_dependent),
1183  "Cannot apply an integration scheme to " << varname);
1184  it->second.ptsc = ptsc;
1185 // if (!first_step)
1186 // GMM_WARNING2("When you apply a scheme to new variable or change the "
1187 // "scheme of a variable after the first time step, "
1188 // "the precomputation of time derivative will not be "
1189 // "executed. Caution: You have to care by yourself of "
1190 // "the compatbility of the operation");
1191  time_integration = 1;
1192  }
1193 
1194  void model::copy_init_time_derivative() {
1195 
1196  for (VAR_SET::iterator it = variables.begin();
1197  it != variables.end(); ++it)
1198  if (it->second.is_variable && it->second.ptsc) {
1199 
1200  std::string name_v, name_previous_v;
1201  it->second.ptsc->time_derivative_to_be_initialized(name_v,
1202  name_previous_v);
1203 
1204  if (name_v.size()) {
1205  if (is_complex()) {
1206  model_complex_plain_vector v0 = complex_variable(name_v);
1207  gmm::copy(v0, set_complex_variable(name_previous_v));
1208  } else {
1209  const model_real_plain_vector &v0 = real_variable(name_v);
1210  gmm::copy(v0, set_real_variable(name_previous_v));
1211  }
1212  }
1213  }
1214  }
1215 
1216  // ----------------------------------------------------------------------
1217  //
1218  // Theta-method scheme for first order problems
1219  //
1220  // ----------------------------------------------------------------------
1221 
1222  class APIDECL first_order_theta_method_scheme
1223  : public virtual_time_scheme {
1224 
1225  std::string U, U0, V, V0;
1226  scalar_type theta;
1227 
1228  public:
1229  // V = (U-U0)/(theta*dt) - ((1-theta)/theta)*V0
1230  virtual void init_affine_dependent_variables(model &md) const {
1231  scalar_type dt = md.get_time_step();
1232  scalar_type a = scalar_type(1)/(theta*dt);
1233  scalar_type b = (scalar_type(1)-theta)/theta;
1234  md.set_factor_of_variable(V, a);
1235  if (md.is_complex()) {
1236  gmm::add(gmm::scaled(md.complex_variable(U0), -complex_type(a)),
1237  gmm::scaled(md.complex_variable(V0), -complex_type(b)),
1238  md.set_complex_constant_part(V));
1239 
1240  } else {
1241  gmm::add(gmm::scaled(md.real_variable(U0), -a),
1242  gmm::scaled(md.real_variable(V0), -b),
1243  md.set_real_constant_part(V));
1244  }
1245  }
1246 
1247  // V = (U-U0)/dt (backward Euler for precomputation)
1248  virtual void init_affine_dependent_variables_precomputation(model &md)
1249  const {
1250  scalar_type dt = md.get_time_step();
1251  md.set_factor_of_variable(V, scalar_type(1)/dt);
1252  if (md.is_complex()) {
1253  gmm::copy(gmm::scaled(md.complex_variable(U0), -complex_type(1)/dt),
1254  md.set_complex_constant_part(V));
1255 
1256  } else {
1257  gmm::copy(gmm::scaled(md.real_variable(U0), -scalar_type(1)/dt),
1258  md.set_real_constant_part(V));
1259  }
1260  }
1261 
1262  virtual void time_derivative_to_be_initialized
1263  (std::string &name_v, std::string &name_previous_v) const
1264  { if (theta != scalar_type(1)) { name_v = V; name_previous_v = V0; } }
1265 
1266  virtual void shift_variables(model &md) const {
1267  if (md.is_complex()) {
1268  gmm::copy(md.complex_variable(U), md.set_complex_variable(U0));
1269  gmm::copy(md.complex_variable(V), md.set_complex_variable(V0));
1270  } else {
1271  gmm::copy(md.real_variable(U), md.set_real_variable(U0));
1272  gmm::copy(md.real_variable(V), md.set_real_variable(V0));
1273  }
1274  }
1275 
1276 
1277  first_order_theta_method_scheme(model &md, std::string varname,
1278  scalar_type th) {
1279  U = varname;
1280  U0 = "Previous_" + U;
1281  V = "Dot_" + U;
1282  V0 = "Previous_Dot_" + U;
1283  theta = th;
1284  GMM_ASSERT1(theta > scalar_type(0) && theta <= scalar_type(1),
1285  "Invalid value of theta parameter for the theta-method");
1286 
1287  if (!(md.variable_exists(V)))
1288  md.add_affine_dependent_variable(V, U);
1289  const mesh_fem *mf = md.pmesh_fem_of_variable(U);
1290  size_type s = md.is_complex() ? gmm::vect_size(md.complex_variable(U))
1291  : gmm::vect_size(md.real_variable(U));
1292 
1293  if (mf) {
1294  if (!(md.variable_exists(U0))) md.add_fem_data(U0, *mf);
1295  if (!(md.variable_exists(V0))) md.add_fem_data(V0, *mf);
1296  } else {
1297  if (!(md.variable_exists(U0))) md.add_fixed_size_data(U0, s);
1298  if (!(md.variable_exists(V0))) md.add_fixed_size_data(V0, s);
1299  }
1300  }
1301 
1302 
1303  };
1304 
1305  void add_theta_method_for_first_order(model &md, const std::string &varname,
1306  scalar_type theta) {
1307  ptime_scheme ptsc
1308  = std::make_shared<first_order_theta_method_scheme>(md, varname,theta);
1309  md.add_time_integration_scheme(varname, ptsc);
1310  }
1311 
1312  // ----------------------------------------------------------------------
1313  //
1314  // Theta-method for second order problems
1315  //
1316  // ----------------------------------------------------------------------
1317 
1318  class APIDECL second_order_theta_method_scheme
1319  : public virtual_time_scheme {
1320 
1321  std::string U, U0, V, V0, A, A0;
1322  scalar_type theta;
1323 
1324  public:
1325  // V = (U-U0)/(theta*dt) - dt*(1-theta)*V0
1326  // A = (U-U0)/(theta^2*dt^2) - V0/(theta^2*dt) - dt*(1-theta)*A0
1327  virtual void init_affine_dependent_variables(model &md) const {
1328  scalar_type dt = md.get_time_step();
1329  md.set_factor_of_variable(V, scalar_type(1)/(theta*dt));
1330  md.set_factor_of_variable(A, scalar_type(1)/(theta*theta*dt*dt));
1331  if (md.is_complex()) {
1332  gmm::add(gmm::scaled(md.complex_variable(U0),
1333  -complex_type(1)/(theta*dt)),
1334  gmm::scaled(md.complex_variable(V0),
1335  -(complex_type(1)-complex_type(theta))/theta),
1336  md.set_complex_constant_part(V));
1337  gmm::add(gmm::scaled(md.complex_variable(U0),
1338  -complex_type(1)/(theta*theta*dt*dt)),
1339  gmm::scaled(md.complex_variable(A0),
1340  -(complex_type(1)-complex_type(theta))/theta),
1341  md.set_complex_constant_part(A));
1342  gmm::add(gmm::scaled(md.complex_variable(V0),
1343  -complex_type(1)/(theta*theta*dt)),
1344  md.set_complex_constant_part(A));
1345 
1346 
1347  } else {
1348  gmm::add(gmm::scaled(md.real_variable(U0),
1349  -scalar_type(1)/(theta*dt)),
1350  gmm::scaled(md.real_variable(V0),
1351  -(scalar_type(1)-theta)/theta),
1352  md.set_real_constant_part(V));
1353  gmm::add(gmm::scaled(md.real_variable(U0),
1354  -scalar_type(1)/(theta*theta*dt*dt)),
1355  gmm::scaled(md.real_variable(A0),
1356  -(scalar_type(1)-theta)/theta),
1357  md.set_real_constant_part(A));
1358  gmm::add(gmm::scaled(md.real_variable(V0),
1359  -scalar_type(1)/(theta*theta*dt)),
1360  md.set_real_constant_part(A));
1361 
1362  }
1363  }
1364 
1365  // V = (U-U0)/dt (backward Euler for precomputation)
1366  // A = (U-U0)/(dt^2) - V0/dt
1367  virtual void init_affine_dependent_variables_precomputation(model &md)
1368  const {
1369  scalar_type dt = md.get_time_step();
1370  md.set_factor_of_variable(V, scalar_type(1)/dt);
1371  md.set_factor_of_variable(A, scalar_type(1)/(dt*dt));
1372  if (md.is_complex()) {
1373  gmm::copy(gmm::scaled(md.complex_variable(U0),
1374  -complex_type(1)/dt),
1375  md.set_complex_constant_part(V));
1376  gmm::add(gmm::scaled(md.complex_variable(U0),
1377  -complex_type(1)/(dt*dt)),
1378  gmm::scaled(md.complex_variable(V0),
1379  -complex_type(1)/dt),
1380  md.set_complex_constant_part(A));
1381  } else {
1382  gmm::copy(gmm::scaled(md.real_variable(U0),
1383  -scalar_type(1)/dt),
1384  md.set_real_constant_part(V));
1385  gmm::add(gmm::scaled(md.real_variable(U0),
1386  -scalar_type(1)/(dt*dt)),
1387  gmm::scaled(md.real_variable(V0),
1388  -scalar_type(1)/dt),
1389  md.set_real_constant_part(A));
1390  }
1391  }
1392 
1393  virtual void time_derivative_to_be_initialized
1394  (std::string &name_v, std::string &name_previous_v) const
1395  { if (theta != scalar_type(1)) { name_v = A; name_previous_v = A0; } }
1396 
1397  virtual void shift_variables(model &md) const {
1398  if (md.is_complex()) {
1399  gmm::copy(md.complex_variable(U), md.set_complex_variable(U0));
1400  gmm::copy(md.complex_variable(V), md.set_complex_variable(V0));
1401  gmm::copy(md.complex_variable(A), md.set_complex_variable(A0));
1402  } else {
1403  gmm::copy(md.real_variable(U), md.set_real_variable(U0));
1404  gmm::copy(md.real_variable(V), md.set_real_variable(V0));
1405  gmm::copy(md.real_variable(A), md.set_real_variable(A0));
1406  }
1407  }
1408 
1409 
1410  second_order_theta_method_scheme(model &md, std::string varname,
1411  scalar_type th) {
1412  U = varname;
1413  U0 = "Previous_" + U;
1414  V = "Dot_" + U;
1415  V0 = "Previous_Dot_" + U;
1416  A = "Dot2_" + U;
1417  A0 = "Previous_Dot2_" + U;
1418  theta = th;
1419  GMM_ASSERT1(theta > scalar_type(0) && theta <= scalar_type(1),
1420  "Invalid value of theta parameter for the theta-method");
1421 
1422  if (!(md.variable_exists(V)))
1423  md.add_affine_dependent_variable(V, U);
1424  if (!(md.variable_exists(A)))
1425  md.add_affine_dependent_variable(A, U);
1426 
1427  const mesh_fem *mf = md.pmesh_fem_of_variable(U);
1428  size_type s = md.is_complex() ? gmm::vect_size(md.complex_variable(U))
1429  : gmm::vect_size(md.real_variable(U));
1430 
1431  if (mf) {
1432  if (!(md.variable_exists(U0))) md.add_fem_data(U0, *mf);
1433  if (!(md.variable_exists(V0))) md.add_fem_data(V0, *mf);
1434  if (!(md.variable_exists(A0))) md.add_fem_data(A0, *mf);
1435  } else {
1436  if (!(md.variable_exists(U0))) md.add_fixed_size_data(U0, s);
1437  if (!(md.variable_exists(V0))) md.add_fixed_size_data(V0, s);
1438  if (!(md.variable_exists(A0))) md.add_fixed_size_data(A0, s);
1439  }
1440  }
1441 
1442 
1443  };
1444 
1445  void add_theta_method_for_second_order(model &md, const std::string &varname,
1446  scalar_type theta) {
1447  ptime_scheme ptsc = std::make_shared<second_order_theta_method_scheme>
1448  (md,varname,theta);
1449  md.add_time_integration_scheme(varname, ptsc);
1450  }
1451 
1452 
1453  // ----------------------------------------------------------------------
1454  //
1455  // Newmark method for second order problems
1456  //
1457  // ----------------------------------------------------------------------
1458 
1459  class APIDECL Newmark_scheme
1460  : public virtual_time_scheme {
1461 
1462  std::string U, U0, V, V0, A, A0;
1463  scalar_type beta, gamma;
1464 
1465  public:
1466  // V = (U-U0)/(theta*dt) - dt*(1-theta)*V0
1467  // A = (U-U0)/(theta^2*dt^2) - V0/(theta^2*dt) - dt*(1-theta)*A0
1468  virtual void init_affine_dependent_variables(model &md) const {
1469  scalar_type dt = md.get_time_step();
1470  scalar_type a0 = scalar_type(1)/(beta*dt*dt), a1 = dt*a0;
1471  scalar_type a2 = (scalar_type(1) - scalar_type(2)*beta)
1472  / (scalar_type(2)*beta);
1473  scalar_type b0 = gamma/(beta*dt), b1 = (beta-gamma)/beta;
1474  scalar_type b2 = dt*(scalar_type(1)-gamma/(scalar_type(2)*beta));
1475 
1476  md.set_factor_of_variable(V, b0);
1477  md.set_factor_of_variable(A, a0);
1478  if (md.is_complex()) {
1479  gmm::add(gmm::scaled(md.complex_variable(U0), -complex_type(b0)),
1480  gmm::scaled(md.complex_variable(V0), complex_type(b1)),
1481  md.set_complex_constant_part(V));
1482  gmm::add(gmm::scaled(md.complex_variable(A0), complex_type(b2)),
1483  md.set_complex_constant_part(V));
1484  gmm::add(gmm::scaled(md.complex_variable(U0), -complex_type(a0)),
1485  gmm::scaled(md.complex_variable(V0), -complex_type(a1)),
1486  md.set_complex_constant_part(A));
1487  gmm::add(gmm::scaled(md.complex_variable(A0), -complex_type(a2)),
1488  md.set_complex_constant_part(A));
1489  } else {
1490  gmm::add(gmm::scaled(md.real_variable(U0), -b0),
1491  gmm::scaled(md.real_variable(V0), b1),
1492  md.set_real_constant_part(V));
1493  gmm::add(gmm::scaled(md.real_variable(A0), b2),
1494  md.set_real_constant_part(V));
1495  gmm::add(gmm::scaled(md.real_variable(U0), -a0),
1496  gmm::scaled(md.real_variable(V0), -a1),
1497  md.set_real_constant_part(A));
1498  gmm::add(gmm::scaled(md.real_variable(A0), -a2),
1499  md.set_real_constant_part(A));
1500 
1501  }
1502  }
1503 
1504  // V = (U-U0)/dt (backward Euler for precomputation)
1505  // A = (U-U0)/(dt^2) - V0/dt
1506  virtual void init_affine_dependent_variables_precomputation(model &md)
1507  const {
1508  scalar_type dt = md.get_time_step();
1509  md.set_factor_of_variable(V, scalar_type(1)/dt);
1510  md.set_factor_of_variable(A, scalar_type(1)/(dt*dt));
1511  if (md.is_complex()) {
1512  gmm::copy(gmm::scaled(md.complex_variable(U0),
1513  -complex_type(1)/dt),
1514  md.set_complex_constant_part(V));
1515  gmm::add(gmm::scaled(md.complex_variable(U0),
1516  -complex_type(1)/(dt*dt)),
1517  gmm::scaled(md.complex_variable(V0),
1518  -complex_type(1)/dt),
1519  md.set_complex_constant_part(A));
1520  } else {
1521  gmm::copy(gmm::scaled(md.real_variable(U0),
1522  -scalar_type(1)/dt),
1523  md.set_real_constant_part(V));
1524  gmm::add(gmm::scaled(md.real_variable(U0),
1525  -scalar_type(1)/(dt*dt)),
1526  gmm::scaled(md.real_variable(V0),
1527  -scalar_type(1)/dt),
1528  md.set_real_constant_part(A));
1529  }
1530  }
1531 
1532  virtual void time_derivative_to_be_initialized
1533  (std::string &name_v, std::string &name_previous_v) const {
1534  if (beta != scalar_type(0.5) || gamma != scalar_type(1))
1535  { name_v = A; name_previous_v = A0; }
1536  }
1537 
1538  virtual void shift_variables(model &md) const {
1539  if (md.is_complex()) {
1540  gmm::copy(md.complex_variable(U), md.set_complex_variable(U0));
1541  gmm::copy(md.complex_variable(V), md.set_complex_variable(V0));
1542  gmm::copy(md.complex_variable(A), md.set_complex_variable(A0));
1543  } else {
1544  gmm::copy(md.real_variable(U), md.set_real_variable(U0));
1545  gmm::copy(md.real_variable(V), md.set_real_variable(V0));
1546  gmm::copy(md.real_variable(A), md.set_real_variable(A0));
1547  }
1548  }
1549 
1550 
1551  Newmark_scheme(model &md, std::string varname,
1552  scalar_type be, scalar_type ga) {
1553  U = varname;
1554  U0 = "Previous_" + U;
1555  V = "Dot_" + U;
1556  V0 = "Previous_Dot_" + U;
1557  A = "Dot2_" + U;
1558  A0 = "Previous_Dot2_" + U;
1559  beta = be; gamma = ga;
1560  GMM_ASSERT1(beta > scalar_type(0) && beta <= scalar_type(1)
1561  && gamma >= scalar_type(0.5) && gamma <= scalar_type(1),
1562  "Invalid parameter values for the Newmark scheme");
1563 
1564  if (!(md.variable_exists(V)))
1565  md.add_affine_dependent_variable(V, U);
1566  if (!(md.variable_exists(A)))
1567  md.add_affine_dependent_variable(A, U);
1568 
1569  const mesh_fem *mf = md.pmesh_fem_of_variable(U);
1570  size_type s = md.is_complex() ? gmm::vect_size(md.complex_variable(U))
1571  : gmm::vect_size(md.real_variable(U));
1572 
1573  if (mf) {
1574  if (!(md.variable_exists(U0))) md.add_fem_data(U0, *mf);
1575  if (!(md.variable_exists(V0))) md.add_fem_data(V0, *mf);
1576  if (!(md.variable_exists(A0))) md.add_fem_data(A0, *mf);
1577  } else {
1578  if (!(md.variable_exists(U0))) md.add_fixed_size_data(U0, s);
1579  if (!(md.variable_exists(V0))) md.add_fixed_size_data(V0, s);
1580  if (!(md.variable_exists(A0))) md.add_fixed_size_data(A0, s);
1581  }
1582  }
1583 
1584 
1585  };
1586 
1587  void add_Newmark_scheme(model &md, const std::string &varname,
1588  scalar_type beta, scalar_type gamma) {
1589  ptime_scheme ptsc = std::make_shared<Newmark_scheme>
1590  (md, varname, beta, gamma);
1591  md.add_time_integration_scheme(varname, ptsc);
1592  }
1593 
1594  // ----------------------------------------------------------------------
1595  //
1596  // Houbolt method
1597  //
1598  // ----------------------------------------------------------------------
1599 
1600  class APIDECL Houbolt_scheme
1601  : public virtual_time_scheme {
1602 
1603  std::string U, U01, U02, U03, V, A;
1604 
1605  public:
1606  // V = 1/(6*dt)*(11*U-18*U01+9*U02-2*U03)
1607  // A = 1/(dt**2)*(2*U-5*U01+4*U02-U03)
1608  virtual void init_affine_dependent_variables(model &md) const {
1609  scalar_type dt = md.get_time_step();
1610  scalar_type a0 = scalar_type(2)/(dt*dt);
1611  scalar_type a1 = scalar_type(5)/(dt*dt);
1612  scalar_type a2 = scalar_type(4)/(dt*dt);
1613  scalar_type a3 = scalar_type(1)/(dt*dt);
1614  scalar_type b0 = scalar_type(11)/(scalar_type(6)*dt);
1615  scalar_type b1 = scalar_type(18)/(scalar_type(6)*dt);
1616  scalar_type b2 = scalar_type(9)/(scalar_type(6)*dt);
1617  scalar_type b3 = scalar_type(2)/(scalar_type(6)*dt);
1618 
1619  md.set_factor_of_variable(V, b0);
1620  md.set_factor_of_variable(A, a0);
1621  if (md.is_complex()) {
1622  gmm::add(gmm::scaled(md.complex_variable(U01), -complex_type(b1)),
1623  gmm::scaled(md.complex_variable(U02), complex_type(b2)),
1624  md.set_complex_constant_part(V));
1625  gmm::add(gmm::scaled(md.complex_variable(U03), -complex_type(b3)),
1626  md.set_complex_constant_part(V));
1627  gmm::add(gmm::scaled(md.complex_variable(U01), -complex_type(a1)),
1628  gmm::scaled(md.complex_variable(U02), complex_type(a2)),
1629  md.set_complex_constant_part(A));
1630  gmm::add(gmm::scaled(md.complex_variable(U03), -complex_type(a3)),
1631  md.set_complex_constant_part(A));
1632  } else {
1633  gmm::add(gmm::scaled(md.real_variable(U01), -b1),
1634  gmm::scaled(md.real_variable(U02), b2),
1635  md.set_real_constant_part(V));
1636  gmm::add(gmm::scaled(md.real_variable(U03), -b3),
1637  md.set_real_constant_part(V));
1638  gmm::add(gmm::scaled(md.real_variable(U01), -a1),
1639  gmm::scaled(md.real_variable(U02), a2),
1640  md.set_real_constant_part(A));
1641  gmm::add(gmm::scaled(md.real_variable(U03), -a3),
1642  md.set_real_constant_part(A));
1643  }
1644  }
1645 
1646  virtual void init_affine_dependent_variables_precomputation(model &md)
1647  const {
1648  (void) md;
1649  }
1650 
1651  virtual void time_derivative_to_be_initialized
1652  (std::string &name_v, std::string &name_previous_v) const {
1653  (void) name_v;
1654  (void) name_previous_v;
1655  }
1656 
1657  virtual void shift_variables(model &md) const {
1658  if (md.is_complex()) {
1659  gmm::copy(md.complex_variable(U02), md.set_complex_variable(U03));
1660  gmm::copy(md.complex_variable(U01), md.set_complex_variable(U02));
1661  gmm::copy(md.complex_variable(U), md.set_complex_variable(U01));
1662  } else {
1663  gmm::copy(md.real_variable(U02), md.set_real_variable(U03));
1664  gmm::copy(md.real_variable(U01), md.set_real_variable(U02));
1665  gmm::copy(md.real_variable(U), md.set_real_variable(U01));
1666  }
1667  }
1668 
1669 
1670  Houbolt_scheme(model &md, std::string varname) {
1671  U = varname;
1672  U01 = "Previous_" + U;
1673  U02 = "Previous2_" + U;
1674  U03 = "Previous3_" + U;
1675  V = "Dot_" + U;
1676  A = "Dot2_" + U;
1677 
1678  if (!(md.variable_exists(V)))
1679  md.add_affine_dependent_variable(V, U);
1680  if (!(md.variable_exists(A)))
1681  md.add_affine_dependent_variable(A, U);
1682 
1683  const mesh_fem *mf = md.pmesh_fem_of_variable(U);
1684  size_type s = md.is_complex() ? gmm::vect_size(md.complex_variable(U))
1685  : gmm::vect_size(md.real_variable(U));
1686 
1687  if (mf) {
1688  if (!(md.variable_exists(U01))) md.add_fem_data(U01, *mf);
1689  if (!(md.variable_exists(U02))) md.add_fem_data(U02, *mf);
1690  if (!(md.variable_exists(U03))) md.add_fem_data(U03, *mf);
1691  } else {
1692  if (!(md.variable_exists(U01))) md.add_fixed_size_data(U01, s);
1693  if (!(md.variable_exists(U02))) md.add_fixed_size_data(U02, s);
1694  if (!(md.variable_exists(U03))) md.add_fixed_size_data(U03, s);
1695  }
1696 
1697  }
1698 
1699  };
1700 
1701  void add_Houbolt_scheme(model &md, const std::string &varname) {
1702  ptime_scheme ptsc = std::make_shared<Houbolt_scheme>
1703  (md, varname);
1704  md.add_time_integration_scheme(varname, ptsc);
1705  }
1706 
1707 
1708 
1709 
1710 
1711 
1712 
1713 
1714 
1715  void model::add_time_dispatcher(size_type ibrick, pdispatcher pdispatch) {
1716  GMM_ASSERT1(valid_bricks[ibrick], "Inexistent brick");
1717  pbrick pbr = bricks[ibrick].pbr;
1718 
1719  bricks[ibrick].pdispatch = pdispatch;
1720 
1721  size_type nbrhs = bricks[ibrick].nbrhs
1722  = std::max(size_type(1), pdispatch->nbrhs());
1723 
1724  gmm::resize(bricks[ibrick].coeffs, nbrhs);
1725 
1726  if (is_complex() && pbr->is_complex()) {
1727  bricks[ibrick].cveclist.resize(nbrhs);
1728  bricks[ibrick].cveclist_sym.resize(nbrhs);
1729  for (size_type k = 1; k < nbrhs; ++k) {
1730  bricks[ibrick].cveclist[k] = bricks[ibrick].cveclist[0];
1731  bricks[ibrick].cveclist_sym[k] = bricks[ibrick].cveclist_sym[0];
1732  }
1733  } else {
1734  bricks[ibrick].rveclist.resize(nbrhs);
1735  bricks[ibrick].rveclist_sym.resize(nbrhs);
1736  for (size_type k = 1; k < nbrhs; ++k) {
1737  bricks[ibrick].rveclist[k] = bricks[ibrick].rveclist[0];
1738  bricks[ibrick].rveclist_sym[k] = bricks[ibrick].rveclist_sym[0];
1739  }
1740  }
1741  }
1742 
1743  const std::string &model::varname_of_brick(size_type ind_brick,
1744  size_type ind_var) {
1745  GMM_ASSERT1(valid_bricks[ind_brick], "Inexistent brick");
1746  GMM_ASSERT1(ind_var < bricks[ind_brick].vlist.size(),
1747  "Inexistent brick variable");
1748  return bricks[ind_brick].vlist[ind_var];
1749  }
1750 
1751  const std::string &model::dataname_of_brick(size_type ind_brick,
1752  size_type ind_data) {
1753  GMM_ASSERT1(valid_bricks[ind_brick], "Inexistent brick");
1754  GMM_ASSERT1(ind_data < bricks[ind_brick].dlist.size(),
1755  "Inexistent brick data");
1756  return bricks[ind_brick].dlist[ind_data];
1757  }
1758 
1759  void model::listbricks(std::ostream &ost, size_type base_id) const {
1760  if (valid_bricks.card() == 0)
1761  ost << "Model with no bricks" << endl;
1762  else {
1763  ost << "List of model bricks:" << endl;
1764  for (dal::bv_visitor i(valid_bricks); !i.finished(); ++i) {
1765  ost << "Brick " << std::setw(3) << std::right << i + base_id
1766  << " " << std::setw(20) << std::right
1767  << bricks[i].pbr->brick_name();
1768  if (!(active_bricks[i])) ost << " (deactivated)";
1769  if (bricks[i].pdispatch) ost << " (dispatched)";
1770  ost << endl << " concerned variables: " << bricks[i].vlist[0];
1771  for (size_type j = 1; j < bricks[i].vlist.size(); ++j)
1772  ost << ", " << bricks[i].vlist[j];
1773  ost << "." << endl;
1774  ost << " brick with " << bricks[i].tlist.size() << " term";
1775  if (bricks[i].tlist.size() > 1) ost << "s";
1776  ost << endl;
1777  // + lister les termes
1778  }
1779  }
1780  }
1781 
1782  // before call to asm_real_tangent_terms or asm_complex_tangent_terms
1783  // from the assembly procedure or a time dispatcher
1784  void model::brick_init(size_type ib, build_version version,
1785  size_type rhs_ind) const {
1786  const brick_description &brick = bricks[ib];
1787  bool cplx = is_complex() && brick.pbr->is_complex();
1788 
1789  // Initialization of vector and matrices.
1790  for (size_type j = 0; j < brick.tlist.size(); ++j) {
1791  const term_description &term = brick.tlist[j];
1792  bool isg = term.is_global;
1793  size_type nbgdof = is_complex() ?
1794  gmm::vect_size(crhs) : gmm::vect_size(rrhs);
1795  size_type nbd1 = isg ? nbgdof : variables[term.var1].size();
1796  size_type nbd2 = isg ? nbgdof : (term.is_matrix_term ?
1797  variables[term.var2].size() : 0);
1798  if (term.is_matrix_term &&
1799  (brick.pbr->is_linear() || (version | BUILD_MATRIX))) {
1800  if (version | BUILD_ON_DATA_CHANGE) {
1801  if (cplx)
1802  gmm::resize(brick.cmatlist[j], nbd1, nbd2);
1803  else
1804  gmm::resize(brick.rmatlist[j], nbd1, nbd2);
1805  } else {
1806  if (cplx)
1807  brick.cmatlist[j] = model_complex_sparse_matrix(nbd1, nbd2);
1808  else
1809  brick.rmatlist[j] = model_real_sparse_matrix(nbd1, nbd2);
1810  }
1811  }
1812  if (brick.pbr->is_linear() || (version | BUILD_RHS)) {
1813  for (size_type k = 0; k < brick.nbrhs; ++k) {
1814  if (cplx) {
1815  if (k == rhs_ind) gmm::clear(brick.cveclist[k][j]);
1816  gmm::resize(brick.cveclist[k][j], nbd1);
1817  if (term.is_symmetric && term.var1.compare(term.var2)) {
1818  if (k == rhs_ind) gmm::clear(brick.cveclist_sym[k][j]);
1819  gmm::resize(brick.cveclist_sym[k][j], nbd2);
1820  }
1821  } else {
1822  if (k == rhs_ind) gmm::clear(brick.rveclist[k][j]);
1823  gmm::resize(brick.rveclist[k][j], nbd1);
1824  if (term.is_symmetric && term.var1.compare(term.var2)) {
1825  if (k == rhs_ind) gmm::clear(brick.rveclist_sym[k][j]);
1826  gmm::resize(brick.rveclist_sym[k][j], nbd2);
1827  }
1828  }
1829  }
1830  }
1831  }
1832  }
1833 
1834  void model::post_to_variables_step(){}
1835 
1836  void model::brick_call(size_type ib, build_version version,
1837  size_type rhs_ind) const
1838  {
1839  const brick_description &brick = bricks[ib];
1840  bool cplx = is_complex() && brick.pbr->is_complex();
1841 
1842  brick_init(ib, version, rhs_ind);
1843 
1844  if (cplx)
1845  {
1846  brick.pbr->complex_pre_assembly_in_serial(*this, ib, brick.vlist,
1847  brick.dlist, brick.mims,
1848  brick.cmatlist,
1849  brick.cveclist[rhs_ind],
1850  brick.cveclist_sym[rhs_ind],
1851  brick.region, version);
1852 
1853  /*distributing the resulting vectors and matrices for individual threads.*/
1854  { //brackets are needed because accumulated_distro has constructor/destructor
1855  //semantics (as in RAII)
1856  accumulated_distro<complex_matlist> cmatlist(brick.cmatlist);
1857  accumulated_distro<complex_veclist> cveclist(brick.cveclist[rhs_ind]);
1858  accumulated_distro<complex_veclist> cveclist_sym(brick.cveclist_sym[rhs_ind]);
1859 
1860  /*running the assembly in parallel*/
1862  brick.pbr->asm_complex_tangent_terms(*this, ib, brick.vlist,
1863  brick.dlist, brick.mims,
1864  cmatlist,
1865  cveclist,
1866  cveclist_sym,
1867  brick.region, version);
1868  )
1869  }
1870  brick.pbr->complex_post_assembly_in_serial(*this, ib, brick.vlist,
1871  brick.dlist, brick.mims,
1872  brick.cmatlist,
1873  brick.cveclist[rhs_ind],
1874  brick.cveclist_sym[rhs_ind],
1875  brick.region, version);
1876 
1877  if (brick.is_update_brick) //contributions of pure update bricks must be deleted
1878  {
1879  for (auto &&mat : brick.cmatlist)
1880  gmm::clear(mat);
1881 
1882  for (auto &&vecs : brick.cveclist)
1883  for (auto &&vec : vecs)
1884  gmm::clear(vec);
1885 
1886  for (auto &&vecs : brick.cveclist_sym)
1887  for (auto &&vec : vecs)
1888  gmm::clear(vec);
1889  }
1890  }
1891  else //not cplx
1892  {
1893  brick.pbr->real_pre_assembly_in_serial(*this, ib, brick.vlist,
1894  brick.dlist, brick.mims,
1895  brick.rmatlist,
1896  brick.rveclist[rhs_ind],
1897  brick.rveclist_sym[rhs_ind],
1898  brick.region, version);
1899  {
1900  /*distributing the resulting vectors and matrices for individual threads.*/
1901  accumulated_distro<real_matlist> rmatlist(brick.rmatlist);
1902  accumulated_distro<real_veclist> rveclist(brick.rveclist[rhs_ind]);
1903  accumulated_distro<real_veclist> rveclist_sym(brick.rveclist_sym[rhs_ind]);
1904 
1905  /*running the assembly in parallel*/
1907  brick.pbr->asm_real_tangent_terms(*this, ib, brick.vlist,
1908  brick.dlist, brick.mims,
1909  rmatlist,
1910  rveclist,
1911  rveclist_sym,
1912  brick.region,
1913  version);
1914  );
1915  }
1916  brick.pbr->real_post_assembly_in_serial(*this, ib, brick.vlist,
1917  brick.dlist, brick.mims,
1918  brick.rmatlist,
1919  brick.rveclist[rhs_ind],
1920  brick.rveclist_sym[rhs_ind],
1921  brick.region, version);
1922 
1923  if (brick.is_update_brick) //contributions of pure update bricks must be deleted
1924  {
1925  for (auto &&mat : brick.rmatlist)
1926  gmm::clear(mat);
1927 
1928  for (auto &&vecs : brick.rveclist)
1929  for (auto &&vec : vecs)
1930  gmm::clear(vec);
1931 
1932  for (auto &&vecs : brick.rveclist_sym)
1933  for (auto &&vec : vecs)
1934  gmm::clear(vec);
1935  }
1936  }
1937  }
1938 
1939 
1940  void model::set_dispatch_coeff() {
1941  for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) {
1942  brick_description &brick = bricks[ib];
1943  if (brick.pdispatch)
1944  brick.pdispatch->set_dispatch_coeff(*this, ib);
1945 
1946  }
1947  }
1948 
1950  context_check(); if (act_size_to_be_done) actualize_sizes();
1951  for (auto && v : variables) v.second.clear_temporaries();
1952 
1953  set_dispatch_coeff();
1954 
1955  for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) {
1956  brick_description &brick = bricks[ib];
1957  if (brick.pdispatch) {
1958  if (is_complex() && brick.pbr->is_complex())
1959  brick.pdispatch->next_complex_iter(*this, ib, brick.vlist,
1960  brick.dlist,
1961  brick.cmatlist, brick.cveclist,
1962  brick.cveclist_sym, true);
1963  else
1964  brick.pdispatch->next_real_iter(*this, ib, brick.vlist, brick.dlist,
1965  brick.rmatlist, brick.rveclist,
1966  brick.rveclist_sym, true);
1967  }
1968  }
1969  }
1970 
1972  context_check(); if (act_size_to_be_done) actualize_sizes();
1973  set_dispatch_coeff();
1974 
1975  for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) {
1976  brick_description &brick = bricks[ib];
1977  if (brick.pdispatch) {
1978  if (is_complex() && brick.pbr->is_complex())
1979  brick.pdispatch->next_complex_iter(*this, ib, brick.vlist,
1980  brick.dlist,
1981  brick.cmatlist, brick.cveclist,
1982  brick.cveclist_sym, false);
1983  else
1984  brick.pdispatch->next_real_iter(*this, ib, brick.vlist, brick.dlist,
1985  brick.rmatlist, brick.rveclist,
1986  brick.rveclist_sym, false);
1987  }
1988  }
1989 
1990  for (auto &&v : variables)
1991  for (size_type i = 1; i < v.second.n_iter; ++i) {
1992  if (is_complex())
1993  gmm::copy(v.second.complex_value[i-1], v.second.complex_value[i]);
1994  else
1995  gmm::copy(v.second.real_value[i-1], v.second.real_value[i]);
1996  v.second.v_num_data[i] = act_counter();
1997  }
1998  }
1999 
2000  bool model::is_var_newer_than_brick(const std::string &varname,
2001  size_type ib, size_type niter) const {
2002  const brick_description &brick = bricks[ib];
2003  var_description &vd = variables[varname];
2004  if (niter == size_type(-1)) niter = vd.default_iter;
2005  return (vd.v_num > brick.v_num || vd.v_num_data[niter] > brick.v_num);
2006  }
2007 
2008  bool model::is_var_mf_newer_than_brick(const std::string &varname,
2009  size_type ib) const {
2010  const brick_description &brick = bricks[ib];
2011  var_description &vd = variables[varname];
2012  return (vd.v_num > brick.v_num);
2013  }
2014 
2015  bool model::is_mim_newer_than_brick(const mesh_im &im,
2016  size_type ib) const {
2017  const brick_description &brick = bricks[ib];
2018  return (im.version_number() > brick.v_num);
2019  }
2020 
2021  void model::define_variable_group(const std::string &group_name,
2022  const std::vector<std::string> &nl) {
2023  GMM_ASSERT1(!(variable_exists(group_name)), "The name of a group of "
2024  "variables cannot be the same as a variable name");
2025 
2026  std::set<const mesh *> ms;
2027  bool is_data_ = false;
2028  for (size_type i = 0; i < nl.size(); ++i) {
2029  if (i == 0)
2030  is_data_ = is_true_data(nl[i]);
2031  else {
2032  GMM_ASSERT1(is_data_ == is_true_data(nl[i]),
2033  "It is not possible to mix variables and data in a group");
2034  }
2035  GMM_ASSERT1(variable_exists(nl[i]),
2036  "All variables in a group have to exist in the model");
2037  const mesh_fem *mf = pmesh_fem_of_variable(nl[i]);
2038  GMM_ASSERT1(mf, "Variables in a group should be fem variables");
2039  GMM_ASSERT1(ms.find(&(mf->linked_mesh())) == ms.end(),
2040  "Two variables in a group cannot share the same mesh");
2041  ms.insert(&(mf->linked_mesh()));
2042  }
2043  variable_groups[group_name] = nl;
2044  }
2045 
2046  void model::add_assembly_assignments(const std::string &varname,
2047  const std::string &expr, size_type rg,
2048  size_type order, bool before) {
2049  GMM_ASSERT1(order < 3 || order == size_type(-1), "Bad order value");
2050  const im_data *imd = pim_data_of_variable(varname);
2051  GMM_ASSERT1(imd != 0, "Only applicable to im_data");
2052  assignement_desc as;
2053  as.varname = varname; as.expr = expr; as.region = rg; as.order = order;
2054  as.before = before;
2055  assignments.push_back(as);
2056  }
2057 
2058  void model::add_temporaries(const varnamelist &vl,
2059  gmm::uint64_type id_num) const {
2060  for (size_type i = 0; i < vl.size(); ++i) {
2061  var_description &vd = variables[vl[i]];
2062  if (vd.n_iter > 1) {
2063  vd.add_temporary(id_num);
2064  }
2065  }
2066  }
2067 
2068  bool model::temporary_uptodate(const std::string &varname,
2069  gmm::uint64_type id_num,
2070  size_type &ind) const {
2071  var_description &vd = variables[varname];
2072  ind = vd.n_iter;
2073  for (; ind < vd.n_iter + vd.n_temp_iter ; ++ind) {
2074  if (vd.v_num_var_iter[ind] == id_num) break;
2075  }
2076  if (ind < vd.n_iter + vd.n_temp_iter) {
2077  if (vd.v_num_iter[ind] <= vd.v_num_data[vd.default_iter]) {
2078  vd.v_num_iter[ind] = act_counter();
2079  return false;
2080  }
2081  return true;
2082  }
2083  ind = size_type(-1);
2084  return true;
2085  }
2086 
2087  void model::set_default_iter_of_variable(const std::string &varname,
2088  size_type ind) const {
2089  if (ind != size_type(-1)) {
2090  var_description &vd = variables[varname];
2091  GMM_ASSERT1(ind < vd.n_iter + vd.n_temp_iter,
2092  "Inexistent iteration " << ind);
2093  vd.default_iter = ind;
2094  }
2095  }
2096 
2097  void model::reset_default_iter_of_variables(const varnamelist &vl) const {
2098  for (size_type i = 0; i < vl.size(); ++i)
2099  variables[vl[i]].default_iter = 0;
2100  }
2101 
2102  const model_real_sparse_matrix &
2103  model::linear_real_matrix_term(size_type ib, size_type iterm) {
2104  GMM_ASSERT1(bricks[ib].tlist[iterm].is_matrix_term,
2105  "Not a matrix term !");
2106  GMM_ASSERT1(bricks[ib].pbr->is_linear(), "Nonlinear term !");
2107  return bricks[ib].rmatlist[iterm];
2108  }
2109 
2110  const model_complex_sparse_matrix &
2111  model::linear_complex_matrix_term(size_type ib, size_type iterm) {
2112  GMM_ASSERT1(bricks[ib].tlist[iterm].is_matrix_term,
2113  "Not a matrix term !");
2114  GMM_ASSERT1(bricks[ib].pbr->is_linear(), "Nonlinear term !");
2115  return bricks[ib].cmatlist[iterm];
2116  }
2117 
2118  // Call the brick to compute the terms
2119  void model::update_brick(size_type ib, build_version version) const {
2120  const brick_description &brick = bricks[ib];
2121  bool cplx = is_complex() && brick.pbr->is_complex();
2122  bool tobecomputed = brick.terms_to_be_computed
2123  || brick.pbr->is_to_be_computed_each_time()
2124  || !(brick.pbr->is_linear());
2125 
2126  // check variable list to test if a mesh_fem has changed.
2127  if (!tobecomputed ) {
2128  for (size_type i = 0; i < brick.vlist.size() && !tobecomputed; ++i) {
2129  var_description &vd = variables[brick.vlist[i]];
2130  if (vd.v_num > brick.v_num) tobecomputed = true;
2131  }
2132  }
2133 
2134  // check data list to test if a vector value of a data has changed.
2135  for (size_type i = 0; i < brick.dlist.size() && !tobecomputed; ++i) {
2136  var_description &vd = variables[brick.dlist[i]];
2137  if (vd.v_num > brick.v_num || vd.v_num_data[vd.default_iter] > brick.v_num) {
2138  tobecomputed = true;
2139  version = build_version(version | BUILD_ON_DATA_CHANGE);
2140  }
2141  }
2142 
2143  // Check if a mesh_im has changed
2144  if (!tobecomputed ) {
2145  for (size_type i = 0; i < brick.mims.size() && !tobecomputed; ++i) {
2146  if (brick.mims[i]->version_number() > brick.v_num) tobecomputed = true;
2147  }
2148  }
2149 
2150  if (tobecomputed) {
2151  brick.external_load = scalar_type(0);
2152 
2153  if (!(brick.pdispatch))
2154  { brick_call(ib, version, 0); }
2155  else {
2156  if (cplx)
2157  brick.pdispatch->asm_complex_tangent_terms
2158  (*this, ib, brick.cmatlist, brick.cveclist, brick.cveclist_sym,
2159  version);
2160  else
2161  brick.pdispatch->asm_real_tangent_terms
2162  (*this, ib, brick.rmatlist, brick.rveclist, brick.rveclist_sym,
2163  version);
2164  }
2165  brick.v_num = act_counter();
2166  }
2167 
2168  if (brick.pbr->is_linear()) brick.terms_to_be_computed = false;
2169  }
2170 
2171  // OBSOLETE (linked to time dispatchers) or to be corrected to take
2172  // into account global matrices
2173  void model::linear_brick_add_to_rhs(size_type ib, size_type ind_data,
2174  size_type n_iter) const {
2175  const brick_description &brick = bricks[ib];
2176  if (brick.pbr->is_linear()) {
2177 
2178  bool cplx = is_complex() && brick.pbr->is_complex();
2179 
2180  for (size_type j = 0; j < brick.tlist.size(); ++j) {
2181  const term_description &term = brick.tlist[j];
2182  bool isg = term.is_global;
2183  size_type nbgdof = nb_dof();
2184 
2185  size_type n_iter_1 = n_iter, n_iter_2 = n_iter;
2186  if (!isg && n_iter == size_type(-1)) {
2187  n_iter_1 = variables[term.var1].default_iter;
2188  if (term.is_matrix_term)
2189  n_iter_2 = variables[term.var2].default_iter;
2190  }
2191 
2192 
2193 
2194  if (term.is_matrix_term) {
2195  if (cplx) {
2196  if (isg) {
2197  model_complex_plain_vector V(nbgdof);
2198  for (VAR_SET::iterator it = variables.begin();
2199  it != variables.end(); ++it)
2200  if (it->second.is_variable) {
2201  size_type n_iter_i = (n_iter == size_type(-1))
2202  ? it->second.default_iter : n_iter;
2203  gmm::copy(it->second.complex_value[n_iter_i],
2204  gmm::sub_vector(V, it->second.I));
2205  }
2207  (brick.cmatlist[j],
2208  gmm::scaled(V, complex_type(-1)),
2209  brick.cveclist[ind_data][j]);
2210  } else
2212  (brick.cmatlist[j],
2213  gmm::scaled(variables[term.var2].complex_value[n_iter_2],
2214  complex_type(-1)),
2215  brick.cveclist[ind_data][j]);
2216  }
2217  else {
2218  if (isg) {
2219  model_real_plain_vector V(nbgdof);
2220  for (VAR_SET::iterator it = variables.begin();
2221  it != variables.end(); ++it)
2222  if (it->second.is_variable) {
2223  size_type n_iter_i = (n_iter == size_type(-1))
2224  ? it->second.default_iter : n_iter;
2225  gmm::copy(it->second.real_value[n_iter_i],
2226  gmm::sub_vector(V, it->second.I));
2227  }
2229  (brick.rmatlist[j], gmm::scaled(V, scalar_type(-1)),
2230  brick.rveclist[ind_data][j]);
2231  } else
2233  (brick.rmatlist[j],
2234  gmm::scaled(variables[term.var2].real_value[n_iter_2],
2235  scalar_type(-1)), brick.rveclist[ind_data][j]);
2236  }
2237 
2238  if (term.is_symmetric && term.var1.compare(term.var2)) {
2239  if (cplx)
2241  (gmm::conjugated(brick.cmatlist[j]),
2242  gmm::scaled(variables[term.var1].complex_value[n_iter_1],
2243  complex_type(-1)),
2244  brick.cveclist_sym[ind_data][j]);
2245  else
2247  (gmm::transposed(brick.rmatlist[j]),
2248  gmm::scaled(variables[term.var1].real_value[n_iter_1],
2249  scalar_type(-1)),
2250  brick.rveclist_sym[ind_data][j]);
2251  }
2252  }
2253  }
2254  }
2255  }
2256 
2257  void model::update_affine_dependent_variables() {
2258  for (VAR_SET::iterator it = variables.begin(); it != variables.end(); ++it)
2259  if (it->second.is_affine_dependent) {
2260  VAR_SET::iterator it2 = variables.find(it->second.org_name);
2261  if (it->second.size() != it2->second.size())
2262  it->second.set_size();
2263  if (it->second.is_complex) {
2264  gmm::add(gmm::scaled(it2->second.complex_value[0],
2265  complex_type(it->second.alpha)),
2266  it->second.affine_complex_value,
2267  it->second.complex_value[0]);
2268  } else {
2269  gmm::add(gmm::scaled(it2->second.real_value[0], it->second.alpha),
2270  it->second.affine_real_value, it->second.real_value[0]);
2271  }
2272  it->second.v_num = std::max(it->second.v_num, it2->second.v_num);
2273  for (size_type i = 0; i < it->second.n_iter; ++i)
2274  {
2275  it->second.v_num_data[i] = std::max(it->second.v_num_data[i],
2276  it2->second.v_num_data[i]);
2277  }
2278  }
2279  }
2280 
2281 
2282 
2283  std::string model::Neumann_term(const std::string &varname,
2284  size_type region) {
2285  std::string result;
2286  const mesh_fem *mf = pmesh_fem_of_variable(varname);
2287  GMM_ASSERT1(mf, "Works only with fem variables.");
2288  mesh &m = const_cast<mesh &>(mf->linked_mesh());
2289  mesh_im dummy_mim(m);
2290 
2291  for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) {
2292  brick_description &brick = bricks[ib];
2293 
2294  bool detected = false;
2295  for (size_type i = 0; i < brick.vlist.size(); ++i)
2296  if (brick.vlist[i].compare(varname) == 0)
2297  { detected = true; break; }
2298 
2299  if (detected && brick.mims.size()) {
2300  int ifo = -1;
2301  for (auto &pmim : brick.mims)
2302  ifo = std::max(ifo, mf->linked_mesh().region(region)
2303  .region_is_faces_of(m, brick.region,
2304  pmim->linked_mesh()));
2305  GMM_ASSERT1(ifo >= 0,
2306  "The given region is only partially covered by "
2307  "region of brick \"" << brick.pbr->brick_name()
2308  << "\". Please subdivise the region");
2309  if (ifo == 1) {
2310  std::string expr = brick.pbr->declare_volume_assembly_string
2311  (*this, ib, brick.vlist, brick.dlist);
2312 
2313  ga_workspace workspace(*this);
2314  size_type order = workspace.add_expression
2315  (expr, dummy_mim, region);
2316  GMM_ASSERT1(order <= 1, "Wrong order for a Neumann term");
2317  expr = workspace.extract_Neumann_term(varname);
2318  if (expr.size()) {
2319  if (result.size())
2320  result += " + " + workspace.extract_Neumann_term(varname);
2321  else
2322  result = workspace.extract_Neumann_term(varname);
2323  }
2324  }
2325  }
2326  }
2327  return result;
2328  }
2329 
2330 
2331 
2332  void model::assembly(build_version version) {
2333 
2334  GMM_ASSERT1(version != BUILD_ON_DATA_CHANGE,
2335  "Invalid assembly version BUILD_ON_DATA_CHANGE");
2336  GMM_ASSERT1(version != BUILD_WITH_LIN,
2337  "Invalid assembly version BUILD_WITH_LIN");
2338  GMM_ASSERT1(version != BUILD_WITH_INTERNAL,
2339  "Invalid assembly version BUILD_WITH_INTERNAL");
2340  int nbp=1;
2341 #if GETFEM_PARA_LEVEL > 1
2342  double t_ref = MPI_Wtime();
2343  int rk=0;
2344  MPI_Comm_rank(MPI_COMM_WORLD, &rk);
2345  MPI_Comm_size(MPI_COMM_WORLD, &nbp);
2346 #endif
2347 
2348  context_check(); if (act_size_to_be_done) actualize_sizes();
2349  if (is_complex()) {
2350  if (version & BUILD_MATRIX) gmm::clear(cTM);
2351  if (version & BUILD_RHS) gmm::clear(crhs);
2352  }
2353  else {
2354  if (version & BUILD_MATRIX) gmm::clear(rTM);
2355  if (version & BUILD_RHS) gmm::clear(rrhs);
2356  }
2357  clear_dof_constraints();
2358  generic_expressions.clear();
2359  update_affine_dependent_variables();
2360 
2361  if (version & BUILD_RHS) approx_external_load_ = scalar_type(0);
2362 
2363  for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) {
2364 
2365  brick_description &brick = bricks[ib];
2366 
2367  // Disables the brick if all its variables are disabled.
2368  bool auto_disabled_brick = true;
2369  for (size_type j = 0; j < brick.vlist.size(); ++j) {
2370  if (!(is_disabled_variable(brick.vlist[j])))
2371  auto_disabled_brick = false;
2372  }
2373  if (auto_disabled_brick) continue;
2374 
2375  update_brick(ib, version);
2376 
2377  bool cplx = is_complex() && brick.pbr->is_complex();
2378 
2379  scalar_type coeff0 = scalar_type(1);
2380  if (brick.pdispatch) coeff0 = brick.matrix_coeff;
2381 
2382  // Assembly of terms
2383 
2384  for (size_type j = 0; j < brick.tlist.size(); ++j) {
2385  term_description &term = brick.tlist[j];
2386  bool isg = term.is_global, isprevious = false;
2387  size_type nbgdof = nb_dof();
2388  scalar_type alpha = coeff0, alpha1 = coeff0, alpha2 = coeff0;
2389  gmm::sub_interval I1(0,nbgdof), I2(0,nbgdof);
2390  var_description *var1=nullptr, *var2=nullptr;
2391  if (!isg) {
2392  VAR_SET::iterator it1 = variables.find(term.var1);
2393  GMM_ASSERT1(it1 != variables.end(), "Internal error");
2394  var1 = &(it1->second);
2395  GMM_ASSERT1(var1->is_variable, "Assembly of data not allowed");
2396  I1 = var1->I;
2397  if (term.is_matrix_term) {
2398  VAR_SET::iterator it2 = variables.find(term.var2);
2399  GMM_ASSERT1(it2 != variables.end(), "Internal error");
2400  var2 = &(it2->second);
2401  I2 = var2->I;
2402  if (!(var2->is_variable)) {
2403  std::string vorgname = sup_previous_and_dot_to_varname(term.var2);
2404  VAR_SET::iterator it3 = variables.find(vorgname);
2405  GMM_ASSERT1(it3->second.is_variable,
2406  "Assembly of data not allowed");
2407  I2 = it3->second.I;
2408  isprevious = true;
2409  }
2410  alpha *= var1->alpha * var2->alpha;
2411  alpha1 *= var1->alpha;
2412  alpha2 *= var2->alpha;
2413  }
2414  }
2415 
2416  if (cplx) { // complex term in complex model
2417  if (term.is_matrix_term && (version & BUILD_MATRIX) && !isprevious
2418  && (isg || (var1->is_enabled() && var2->is_enabled()))) {
2419  gmm::add(gmm::scaled(brick.cmatlist[j], alpha),
2420  gmm::sub_matrix(cTM, I1, I2));
2421  if (term.is_symmetric && I1.first() != I2.first())
2422  gmm::add(gmm::scaled(gmm::transposed(brick.cmatlist[j]), alpha),
2423  gmm::sub_matrix(cTM, I2, I1));
2424  }
2425  if (version & BUILD_RHS) {
2426  //FIXME MPI_SUM_VECTOR(crhs)
2427  if (isg || var1->is_enabled()) {
2428  if (brick.pdispatch)
2429  for (size_type k = 0; k < brick.nbrhs; ++k)
2430  gmm::add(gmm::scaled(brick.cveclist[k][j],
2431  brick.coeffs[k]),
2432  gmm::sub_vector(crhs, I1));
2433  else
2434  gmm::add(gmm::scaled(brick.cveclist[0][j],
2435  complex_type(alpha1)),
2436  gmm::sub_vector(crhs, I1));
2437  }
2438  if (term.is_matrix_term && brick.pbr->is_linear() && is_linear()) {
2439  if (var2->is_affine_dependent && var1->is_enabled())
2440  gmm::mult_add(brick.cmatlist[j],
2441  gmm::scaled(var2->affine_complex_value,
2442  complex_type(-alpha1)),
2443  gmm::sub_vector(crhs, I1));
2444  if (term.is_symmetric && I1.first() != I2.first()
2445  && var1->is_affine_dependent && var2->is_enabled())
2446  gmm::mult_add(gmm::conjugated(brick.cmatlist[j]),
2447  gmm::scaled(var1->affine_complex_value,
2448  complex_type(-alpha2)),
2449  gmm::sub_vector(crhs, I2));
2450  }
2451  if (term.is_matrix_term && brick.pbr->is_linear()
2452  && (!is_linear() || (version & BUILD_WITH_LIN))
2453  && var1->is_enabled())
2454  gmm::mult_add(brick.cmatlist[j],
2455  gmm::scaled(var2->complex_value[0],
2456  complex_type(-alpha1)),
2457  gmm::sub_vector(crhs, I1));
2458  if (term.is_symmetric && I1.first() != I2.first()
2459  && var2->is_enabled()) {
2460  if (brick.pdispatch)
2461  for (size_type k = 0; k < brick.nbrhs; ++k)
2462  gmm::add(gmm::scaled(brick.cveclist_sym[k][j],
2463  brick.coeffs[k]),
2464  gmm::sub_vector(crhs, I2));
2465  else
2466  gmm::add(gmm::scaled(brick.cveclist_sym[0][j],
2467  complex_type(alpha2)),
2468  gmm::sub_vector(crhs, I2));
2469  if (brick.pbr->is_linear()
2470  && (!is_linear() || (version & BUILD_WITH_LIN)))
2471  gmm::mult_add(gmm::conjugated(brick.cmatlist[j]),
2472  gmm::scaled(var1->complex_value[0],
2473  complex_type(-alpha2)),
2474  gmm::sub_vector(crhs, I2));
2475  }
2476  }
2477  } else if (is_complex()) { // real term in complex model
2478  if (term.is_matrix_term && (version & BUILD_MATRIX) && !isprevious
2479  && (isg || (var1->is_enabled() && var2->is_enabled()))) {
2480  gmm::add(gmm::scaled(brick.rmatlist[j], alpha),
2481  gmm::sub_matrix(cTM, I1, I2));
2482  if (term.is_symmetric && I1.first() != I2.first())
2483  gmm::add(gmm::scaled(gmm::transposed(brick.rmatlist[j]), alpha),
2484  gmm::sub_matrix(cTM, I2, I1));
2485  }
2486  if (version & BUILD_RHS) {
2487  //FIXME MPI_SUM_VECTOR(crhs)
2488  if (isg || var1->is_enabled()) {
2489  if (brick.pdispatch)
2490  for (size_type k = 0; k < brick.nbrhs; ++k)
2491  gmm::add(gmm::scaled(brick.rveclist[k][j],
2492  brick.coeffs[k]),
2493  gmm::sub_vector(crhs, I1));
2494  else
2495  gmm::add(gmm::scaled(brick.rveclist[0][j], alpha1),
2496  gmm::sub_vector(crhs, I1));
2497  }
2498  if (term.is_matrix_term && brick.pbr->is_linear() && is_linear()) {
2499  if (var2->is_affine_dependent && var1->is_enabled())
2500  gmm::mult_add(brick.rmatlist[j],
2501  gmm::scaled(var2->affine_complex_value,
2502  complex_type(-alpha1)),
2503  gmm::sub_vector(crhs, I1));
2504  if (term.is_symmetric && I1.first() != I2.first()
2505  && var1->is_affine_dependent && var2->is_enabled())
2506  gmm::mult_add(gmm::transposed(brick.rmatlist[j]),
2507  gmm::scaled(var1->affine_complex_value,
2508  complex_type(-alpha2)),
2509  gmm::sub_vector(crhs, I2));
2510  }
2511  if (term.is_matrix_term && brick.pbr->is_linear()
2512  && (!is_linear() || (version & BUILD_WITH_LIN))
2513  && var1->is_enabled())
2514  gmm::mult_add(brick.rmatlist[j],
2515  gmm::scaled(var2->complex_value[0],
2516  complex_type(-alpha1)),
2517  gmm::sub_vector(crhs, I1));
2518  if (term.is_symmetric && I1.first() != I2.first()
2519  && var2->is_enabled()) {
2520  if (brick.pdispatch)
2521  for (size_type k = 0; k < brick.nbrhs; ++k)
2522  gmm::add(gmm::scaled(brick.rveclist_sym[k][j],
2523  brick.coeffs[k]),
2524  gmm::sub_vector(crhs, I2));
2525  else
2526  gmm::add(gmm::scaled(brick.rveclist_sym[0][j], alpha2),
2527  gmm::sub_vector(crhs, I2));
2528 
2529  if (brick.pbr->is_linear()
2530  && (!is_linear() || (version & BUILD_WITH_LIN)))
2531  gmm::mult_add(gmm::transposed(brick.rmatlist[j]),
2532  gmm::scaled(var1->complex_value[0],
2533  complex_type(-alpha2)),
2534  gmm::sub_vector(crhs, I2));
2535  }
2536  }
2537  } else { // real term in real model
2538  if (term.is_matrix_term && (version & BUILD_MATRIX) && !isprevious
2539  && (isg || (var1->is_enabled() && var2->is_enabled()))) {
2540  gmm::add(gmm::scaled(brick.rmatlist[j], alpha),
2541  gmm::sub_matrix(rTM, I1, I2));
2542  if (term.is_symmetric && I1.first() != I2.first())
2543  gmm::add(gmm::scaled(gmm::transposed(brick.rmatlist[j]), alpha),
2544  gmm::sub_matrix(rTM, I2, I1));
2545  }
2546  if (version & BUILD_RHS) {
2547  // Contributions to interval I1 of var1
2548  auto vec_out1 = gmm::sub_vector(rrhs, I1);
2549  if (isg || var1->is_enabled()) {
2550  if (brick.pdispatch)
2551  for (size_type k = 0; k < brick.nbrhs; ++k)
2552  gmm::add(gmm::scaled(brick.rveclist[k][j],
2553  brick.coeffs[k]),
2554  vec_out1);
2555  else
2556  gmm::add(gmm::scaled(brick.rveclist[0][j], alpha1),
2557  vec_out1);
2558  }
2559  if (var1->is_enabled()
2560  && term.is_matrix_term && brick.pbr->is_linear()) {
2561  bool affine_contrib(is_linear() && var2->is_affine_dependent);
2562  bool linear_contrib(!is_linear() || (version & BUILD_WITH_LIN));
2563  const auto &matj = brick.rmatlist[j];
2564  const auto vec_affine2 = gmm::scaled(var2->affine_real_value,
2565  -alpha1);
2566  const auto vec_linear2 = gmm::scaled(var2->real_value[0],
2567  -alpha1);
2568  if (nbp > 1) {
2569  model_real_plain_vector vec_tmp1(I1.size(), 0.);
2570  if (affine_contrib) // Affine dependent variable contribution
2571  gmm::mult(matj, vec_affine2, vec_tmp1);
2572  if (linear_contrib) // Linear term contribution
2573  gmm::mult_add(matj, vec_linear2, vec_tmp1);
2574  MPI_SUM_VECTOR(vec_tmp1);
2575  gmm::add(vec_tmp1, vec_out1);
2576  } else { // nbp == 1
2577  if (affine_contrib) // Affine dependent variable contribution
2578  gmm::mult_add(matj, vec_affine2, vec_out1);
2579  if (linear_contrib) // Linear term contribution
2580  gmm::mult_add(matj, vec_linear2, vec_out1);
2581  }
2582  }
2583 
2584  // Contributions to interval I2 of var2 due to symmetric terms
2585  if (term.is_symmetric && I1.first() != I2.first() &&
2586  var2->is_enabled()) {
2587  auto vec_out2 = gmm::sub_vector(rrhs, I2);
2588  if (brick.pdispatch)
2589  for (size_type k = 0; k < brick.nbrhs; ++k)
2590  gmm::add(gmm::scaled(brick.rveclist_sym[k][j],
2591  brick.coeffs[k]),
2592  vec_out2);
2593  else
2594  gmm::add(gmm::scaled(brick.rveclist_sym[0][j], alpha2),
2595  vec_out2);
2596  if (term.is_matrix_term && brick.pbr->is_linear()) {
2597  bool affine_contrib(is_linear() && var1->is_affine_dependent);
2598  bool linear_contrib(!is_linear() || (version & BUILD_WITH_LIN));
2599  const auto matj_trans = gmm::transposed(brick.rmatlist[j]);
2600  const auto vec_affine1 = gmm::scaled(var1->affine_real_value,
2601  -alpha2);
2602  const auto vec_linear1 = gmm::scaled(var1->real_value[0],
2603  -alpha2);
2604  if (nbp > 1) {
2605  model_real_plain_vector vec_tmp2(I2.size(),0.);
2606  if (affine_contrib) // Affine dependent variable contribution
2607  gmm::mult(matj_trans, vec_affine1, vec_tmp2);
2608  if (linear_contrib) // Linear term contribution
2609  gmm::mult_add(matj_trans, vec_linear1, vec_tmp2);
2610  MPI_SUM_VECTOR(vec_tmp2);
2611  gmm::add(vec_tmp2, vec_out2);
2612  } else { // nbp == 1
2613  if (affine_contrib) // Affine dependent variable contribution
2614  gmm::mult_add(matj_trans, vec_affine1, vec_out2);
2615  if (linear_contrib) // Linear term contribution
2616  gmm::mult_add(matj_trans, vec_linear1, vec_out2);
2617  }
2618  }
2619  }
2620  }
2621  }
2622  }
2623 
2624  if (brick.pbr->is_linear())
2625  brick.terms_to_be_computed = false;
2626  // Commented to allow to get the information after assembly. Used in
2627  // some aplications. Should be optional ?
2628 // else
2629 // if (cplx) {
2630 // brick.cmatlist = complex_matlist(brick.tlist.size());
2631 // brick.cveclist[0] = complex_veclist(brick.tlist.size());
2632 // } else {
2633 // brick.rmatlist = real_matlist(brick.tlist.size());
2634 // brick.rveclist[0] = real_veclist(brick.tlist.size());
2635 // }
2636 
2637  if (version & BUILD_RHS) approx_external_load_ += brick.external_load;
2638  }
2639 
2640  if (version & BUILD_RHS && version & BUILD_WITH_INTERNAL) {
2641  GMM_ASSERT1(gmm::vect_size(full_rrhs) > 0 && has_internal_variables(),
2642  "Internal error");
2643  gmm::sub_interval IP(0,gmm::vect_size(rrhs));
2644  gmm::fill(full_rrhs, 0.);
2645  gmm::copy(rrhs, gmm::sub_vector(full_rrhs, IP)); // TICTIC
2646  }
2647 
2648  // Generic expressions
2649  if (generic_expressions.size()) {
2650  GMM_ASSERT1(!is_complex(), "to be done");
2651 
2652  if (version & BUILD_RHS)
2653  GMM_TRACE2("Global generic assembly RHS");
2654  if (version & BUILD_MATRIX)
2655  GMM_TRACE2("Global generic assembly tangent term");
2656 
2657  // auxilliary lambda function
2658  auto add_assignments_and_expressions_to_workspace =
2659  [&](ga_workspace &workspace)
2660  {
2661  for (const auto &ad : assignments)
2662  workspace.add_assignment_expression
2663  (ad.varname, ad.expr, ad.region, ad.order, ad.before);
2664  for (const auto &ge : generic_expressions)
2665  workspace.add_expression(ge.expr, ge.mim, ge.region,
2666  2, ge.secondary_domain);
2667  };
2668 
2669  const bool with_internal = version & BUILD_WITH_INTERNAL
2671  model_real_sparse_matrix intern_mat; // temp for extracting condensation info
2672  model_real_plain_vector res0, // holds the original RHS
2673  res1; // holds the condensed RHS
2674 
2675  size_type full_size = gmm::vect_size(full_rrhs),
2676  primary_size = gmm::vect_size(rrhs);
2677 
2678  if ((version & BUILD_RHS) || (version & BUILD_MATRIX && with_internal))
2679  gmm::resize(res0, with_internal ? full_size : primary_size);
2680  if (version & BUILD_MATRIX && with_internal)
2681  gmm::resize(res1, full_size);
2682 
2683  if (version & BUILD_MATRIX) {
2684  if (with_internal) {
2685  gmm::resize(intern_mat, full_size, primary_size);
2686  gmm::resize(res1, full_size);
2687  }
2688  accumulated_distro<decltype(rTM)> tangent_matrix_distro(rTM);
2689  accumulated_distro<decltype(intern_mat)> intern_mat_distro(intern_mat);
2691 
2692  if (version & BUILD_RHS) { // both BUILD_RHS & BUILD_MATRIX
2694  GETFEM_OMP_PARALLEL( // running the assembly in parallel
2695  ga_workspace workspace(*this);
2696  add_assignments_and_expressions_to_workspace(workspace);
2697  workspace.set_assembled_vector(res0_distro);
2698  workspace.assembly(1, with_internal);
2699  if (with_internal) { // Condensation reads from/writes to rhs
2700  gmm::copy(res0_distro.get(), res1_distro.get());
2701  gmm::add(gmm::scaled(full_rrhs, scalar_type(-1)),
2702  res1_distro.get()); // initial value residual=-rhs (actually only the internal variables residual is needed)
2703  workspace.set_assembled_vector(res1_distro);
2704  workspace.set_internal_coupling_matrix(intern_mat_distro);
2705  }
2706  workspace.set_assembled_matrix(tangent_matrix_distro);
2707  workspace.assembly(2, with_internal);
2708  ) // end GETFEM_OMP_PARALLEL
2709  } // end of res0_distro scope
2710  else { // only BUILD_MATRIX
2711  GETFEM_OMP_PARALLEL( // running the assembly in parallel
2712  ga_workspace workspace(*this);
2713  add_assignments_and_expressions_to_workspace(workspace);
2714  if (with_internal) { // Condensation reads from/writes to rhs
2715  gmm::copy(gmm::scaled(full_rrhs, scalar_type(-1)),
2716  res1_distro.get()); // initial value residual=-rhs (actually only the internal variables residual is needed)
2717  workspace.set_assembled_vector(res1_distro);
2718  workspace.set_internal_coupling_matrix(intern_mat_distro);
2719  }
2720  workspace.set_assembled_matrix(tangent_matrix_distro);
2721  workspace.assembly(2, with_internal);
2722  ) // end GETFEM_OMP_PARALLEL
2723  }
2724  } // end of tangent_matrix_distro, intern_mat_distro, res1_distro scope
2725  else if (version & BUILD_RHS) {
2727  GETFEM_OMP_PARALLEL( // running the assembly in parallel
2728  ga_workspace workspace(*this);
2729  add_assignments_and_expressions_to_workspace(workspace);
2730  workspace.set_assembled_vector(res0_distro);
2731  workspace.assembly(1, with_internal);
2732  ) // end GETFEM_OMP_PARALLEL
2733  } // end of res0_distro scope
2734 
2735  if (version & BUILD_RHS) {
2736  gmm::scale(res0, scalar_type(-1)); // from residual to rhs
2737  if (with_internal) {
2738  gmm::sub_interval IP(0,gmm::vect_size(rrhs));
2739  gmm::add(gmm::sub_vector(res0, IP), rrhs); // TOCTOC
2740  gmm::add(res0, full_rrhs);
2741  } else
2742  gmm::add(res0, rrhs);
2743  }
2744 
2745  if (version & BUILD_MATRIX && with_internal) {
2746  gmm::scale(res1, scalar_type(-1)); // from residual to rhs
2747  gmm::sub_interval IP(0, primary_size),
2748  II(primary_size, full_size-primary_size);
2749  gmm::copy(gmm::sub_matrix(intern_mat, II, IP), internal_rTM); // --> internal_rTM
2750  gmm::add(gmm::sub_vector(res1, IP), rrhs); // --> rrhs
2751  gmm::copy(gmm::sub_vector(res1, II), internal_sol); // --> internal_sol
2752  }
2753  }
2754 
2755  // Post simplification for dof constraints
2756  if ((version & BUILD_RHS) || (version & BUILD_MATRIX)) {
2757  if (is_complex()) {
2758  std::vector<size_type> dof_indices;
2759  std::vector<complex_type> dof_pr_values;
2760  std::vector<complex_type> dof_go_values;
2761  for (const auto &keyval : complex_dof_constraints) {
2762  const gmm::sub_interval &I = interval_of_variable(keyval.first);
2763  const model_complex_plain_vector &V = complex_variable(keyval.first);
2764  for (const auto &val : keyval.second) {
2765  dof_indices.push_back(val.first + I.first());
2766  dof_go_values.push_back(val.second);
2767  dof_pr_values.push_back(V[val.first]);
2768  }
2769  }
2770 
2771  if (dof_indices.size()) {
2772  gmm::sub_index SI(dof_indices);
2773  gmm::sub_interval II(0, nb_dof());
2774 
2775  if (version & BUILD_RHS) {
2776  if (MPI_IS_MASTER())
2777  approx_external_load_ += gmm::vect_norm1(dof_go_values);
2778  if (is_linear_) {
2779  if (is_symmetric_) {
2780  scalar_type valnorm = gmm::vect_norm2(dof_go_values);
2781  if (valnorm > scalar_type(0)) {
2782  GMM_ASSERT1(version & BUILD_MATRIX, "Rhs only for a "
2783  "symmetric linear problem with dof "
2784  "constraint not allowed");
2785  model_complex_plain_vector vv(gmm::vect_size(crhs));
2786  gmm::mult(gmm::sub_matrix(cTM, II, SI), dof_go_values, vv);
2787  MPI_SUM_VECTOR(vv);
2788  gmm::add(gmm::scaled(vv, scalar_type(-1)), crhs);
2789  }
2790  }
2791  gmm::copy(dof_go_values, gmm::sub_vector(crhs, SI));
2792  } else {
2793  gmm::add(dof_go_values,
2794  gmm::scaled(dof_pr_values, complex_type(-1)),
2795  gmm::sub_vector(crhs, SI));
2796  }
2797  }
2798  if (version & BUILD_MATRIX) {
2799  gmm::clear(gmm::sub_matrix(cTM, SI, II));
2800  if (is_symmetric_) gmm::clear(gmm::sub_matrix(cTM, II, SI));
2801 
2802  if (MPI_IS_MASTER()) {
2803  for (size_type i = 0; i < dof_indices.size(); ++i)
2804  cTM(dof_indices[i], dof_indices[i]) = complex_type(1);
2805  }
2806  }
2807  }
2808  } else { // !is_complex()
2809  std::vector<size_type> dof_indices;
2810  std::vector<scalar_type> dof_pr_values;
2811  std::vector<scalar_type> dof_go_values;
2812  for (const auto &keyval : real_dof_constraints) {
2813  const gmm::sub_interval &I = interval_of_variable(keyval.first);
2814  const model_real_plain_vector &V = real_variable(keyval.first);
2815  for (const auto &val : keyval.second) {
2816  dof_indices.push_back(val.first + I.first());
2817  dof_go_values.push_back(val.second);
2818  dof_pr_values.push_back(V[val.first]);
2819  }
2820  }
2821 
2822  #if GETFEM_PARA_LEVEL > 1
2823  GMM_ASSERT1(MPI_IS_MASTER() || (dof_indices.size() == 0),
2824  "Sorry, for the moment, the dof constraints have to be "
2825  "added on the master process only");
2826  size_type dof_indices_size = dof_indices.size();
2827  MPI_BCAST0_SCALAR(dof_indices_size);
2828  dof_indices.resize(dof_indices_size);
2829  MPI_BCAST0_VECTOR(dof_indices);
2830  dof_pr_values.resize(dof_indices_size);
2831  MPI_BCAST0_VECTOR(dof_pr_values);
2832  dof_go_values.resize(dof_indices_size);
2833  MPI_BCAST0_VECTOR(dof_go_values);
2834  #endif
2835 
2836  if (dof_indices.size()) {
2837  gmm::sub_index SI(dof_indices);
2838  gmm::sub_interval II(0, nb_dof());
2839 
2840  if (version & BUILD_RHS) {
2841  if (MPI_IS_MASTER())
2842  approx_external_load_ += gmm::vect_norm1(dof_go_values);
2843  if (is_linear_) {
2844  if (is_symmetric_) {
2845  scalar_type valnorm = gmm::vect_norm2(dof_go_values);
2846  if (valnorm > scalar_type(0)) {
2847  GMM_ASSERT1(version & BUILD_MATRIX, "Rhs only for a "
2848  "symmetric linear problem with dof "
2849  "constraint not allowed");
2850  model_real_plain_vector vv(gmm::vect_size(rrhs));
2851  gmm::mult(gmm::sub_matrix(rTM, II, SI), dof_go_values, vv);
2852  MPI_SUM_VECTOR(vv);
2853  gmm::add(gmm::scaled(vv, scalar_type(-1)), rrhs);
2854  }
2855  }
2856  gmm::copy(dof_go_values, gmm::sub_vector(rrhs, SI));
2857  } else {
2858  gmm::add(dof_go_values,
2859  gmm::scaled(dof_pr_values, scalar_type(-1)),
2860  gmm::sub_vector(rrhs, SI));
2861  }
2862  }
2863  if (version & BUILD_MATRIX) {
2864  gmm::clear(gmm::sub_matrix(rTM, SI, II));
2865  if (is_symmetric_) gmm::clear(gmm::sub_matrix(rTM, II, SI));
2866 
2867  if (MPI_IS_MASTER()) {
2868  for (size_type i = 0; i < dof_indices.size(); ++i)
2869  rTM(dof_indices[i], dof_indices[i]) = scalar_type(1);
2870  }
2871  }
2872  }
2873  }
2874  }
2875 
2876  if (version & BUILD_RHS) {
2877  // some contributions are added only in the master process
2878  // send the correct result to all other processes
2879  MPI_BCAST0_SCALAR(approx_external_load_);
2880  }
2881 
2882  #if GETFEM_PARA_LEVEL > 1
2883  // int rk; MPI_Comm_rank(MPI_COMM_WORLD, &rk);
2884  if (MPI_IS_MASTER()) cout << "Assembly time " << MPI_Wtime()-t_ref << endl;
2885  #endif
2886 
2887  }
2888 
2889 
2890  const mesh_fem &
2891  model::mesh_fem_of_variable(const std::string &name) const {
2892  auto it = find_variable(no_old_prefix_name(name));
2893  return it->second.associated_mf();
2894  }
2895 
2896  const mesh_fem *
2897  model::pmesh_fem_of_variable(const std::string &name) const {
2898  auto it = find_variable(no_old_prefix_name(name));
2899  return it->second.passociated_mf();
2900  }
2901 
2902  bgeot::multi_index
2903  model::qdims_of_variable(const std::string &name) const {
2904  auto it = find_variable(no_old_prefix_name(name));
2905  const mesh_fem *mf = it->second.passociated_mf();
2906  const im_data *imd = it->second.imd;
2907  size_type n = it->second.qdim();
2908  if (mf) {
2909  bgeot::multi_index mi = mf->get_qdims();
2910  if (n > 1 || it->second.qdims.size() > 1) {
2911  size_type i = 0;
2912  if (mi.back() == 1) { mi.back() *= it->second.qdims[0]; ++i; }
2913  for (; i < it->second.qdims.size(); ++i)
2914  mi.push_back(it->second.qdims[i]);
2915  }
2916  return mi;
2917  } else if (imd) {
2918  bgeot::multi_index mi = imd->tensor_size();
2919  if (n > 1 || it->second.qdims.size() > 1) {
2920  size_type i = 0;
2921  if (mi.back() == 1) { mi.back() *= it->second.qdims[0]; ++i; }
2922  for (; i < it->second.qdims.size(); ++i)
2923  mi.push_back(it->second.qdims[i]);
2924  }
2925  return mi;
2926  }
2927  return it->second.qdims;
2928  }
2929 
2930  size_type model::qdim_of_variable(const std::string &name) const {
2931  auto it = find_variable(no_old_prefix_name(name));
2932  const mesh_fem *mf = it->second.passociated_mf();
2933  const im_data *imd = it->second.imd;
2934  size_type n = it->second.qdim();
2935  if (mf) {
2936  return mf->get_qdim() * n;
2937  } else if (imd) {
2938  return imd->tensor_size().total_size() * n;
2939  }
2940  return n;
2941  }
2942 
2943 
2944  const gmm::sub_interval &
2945  model::interval_of_variable(const std::string &name) const {
2946  context_check();
2947  if (act_size_to_be_done) actualize_sizes();
2948  VAR_SET::const_iterator it = find_variable(name);
2949  return it->second.I;
2950  }
2951 
2952  const model_real_plain_vector &
2953  model::real_variable(const std::string &name) const {
2954  return is_old(name) ? real_variable(no_old_prefix_name(name), 1)
2955  : real_variable(name, size_type(-1));
2956  }
2957 
2958  const model_real_plain_vector &
2959  model::real_variable(const std::string &name, size_type niter) const {
2960  GMM_ASSERT1(!complex_version, "This model is a complex one");
2961  GMM_ASSERT1(!is_old(name), "Please don't use Old_ prefix in combination "
2962  "with variable version");
2963  context_check();
2964  auto it = variables.find(name);
2965  GMM_ASSERT1(it != variables.end(), "Undefined variable " << name);
2966  if (act_size_to_be_done && it->second.is_fem_dofs) {
2967  if (it->second.filter != VDESCRFILTER_NO)
2968  actualize_sizes();
2969  else
2970  it->second.set_size();
2971  }
2972  if (niter == size_type(-1)) niter = it->second.default_iter;
2973  GMM_ASSERT1(it->second.n_iter + it->second.n_temp_iter > niter,
2974  "Invalid iteration number " << niter << " for " << name);
2975  return it->second.real_value[niter];
2976  }
2977 
2978  const model_complex_plain_vector &
2979  model::complex_variable(const std::string &name) const {
2980  return is_old(name) ? complex_variable(no_old_prefix_name(name), 1)
2981  : complex_variable(name, size_type(-1));
2982  }
2983 
2984  const model_complex_plain_vector &
2985  model::complex_variable(const std::string &name, size_type niter) const {
2986  GMM_ASSERT1(complex_version, "This model is a real one");
2987  GMM_ASSERT1(!is_old(name), "Please don't use Old_ prefix in combination with"
2988  " variable version");
2989  context_check();
2990  auto it = variables.find(name);
2991  GMM_ASSERT1(it!=variables.end(), "Undefined variable " << name);
2992  if (act_size_to_be_done && it->second.is_fem_dofs) {
2993  if (it->second.filter != VDESCRFILTER_NO)
2994  actualize_sizes();
2995  else
2996  it->second.set_size();
2997  }
2998  if (niter == size_type(-1)) niter = it->second.default_iter;
2999  GMM_ASSERT1(it->second.n_iter + it->second.n_temp_iter > niter,
3000  "Invalid iteration number "
3001  << niter << " for " << name);
3002  return it->second.complex_value[niter];
3003  }
3004 
3005  model_real_plain_vector &
3006  model::set_real_variable(const std::string &name) const {
3007  return is_old(name) ? set_real_variable(no_old_prefix_name(name), 1)
3008  : set_real_variable(name, size_type(-1));
3009  }
3010 
3011 
3012  model_real_plain_vector &
3013  model::set_real_variable(const std::string &name, size_type niter) const {
3014  GMM_ASSERT1(!complex_version, "This model is a complex one");
3015  GMM_ASSERT1(!is_old(name), "Please don't use Old_ prefix in combination with"
3016  " variable version");
3017  context_check();
3018  auto it = variables.find(name);
3019  GMM_ASSERT1(it!=variables.end(), "Undefined variable " << name);
3020  if (act_size_to_be_done && it->second.is_fem_dofs) {
3021  if (it->second.filter != VDESCRFILTER_NO)
3022  actualize_sizes();
3023  else
3024  it->second.set_size();
3025  }
3026  if (niter == size_type(-1)) niter = it->second.default_iter;
3027  it->second.v_num_data[niter] = act_counter();
3028  GMM_ASSERT1(it->second.n_iter + it->second.n_temp_iter > niter,
3029  "Invalid iteration number "
3030  << niter << " for " << name);
3031  return it->second.real_value[niter];
3032  }
3033 
3034  model_complex_plain_vector &
3035  model::set_complex_variable(const std::string &name) const {
3036  return is_old(name) ? set_complex_variable(no_old_prefix_name(name), 1)
3037  : set_complex_variable(name, size_type(-1));
3038  }
3039 
3040  model_complex_plain_vector &
3041  model::set_complex_variable(const std::string &name, size_type niter) const {
3042  GMM_ASSERT1(complex_version, "This model is a real one");
3043  GMM_ASSERT1(!is_old(name), "Please don't use Old_ prefix in combination with"
3044  " variable version");
3045  context_check();
3046  auto it = variables.find(name);
3047  GMM_ASSERT1(it!=variables.end(), "Undefined variable " << name);
3048  if (act_size_to_be_done && it->second.is_fem_dofs) {
3049  if (it->second.filter != VDESCRFILTER_NO)
3050  actualize_sizes();
3051  else
3052  it->second.set_size();
3053  }
3054  if (niter == size_type(-1)) niter = it->second.default_iter;
3055  it->second.v_num_data[niter] = act_counter();
3056  GMM_ASSERT1(it->second.n_iter + it->second.n_temp_iter > niter,
3057  "Invalid iteration number "
3058  << niter << " for " << name);
3059  return it->second.complex_value[niter];
3060  }
3061 
3062  model_real_plain_vector &
3063  model::set_real_constant_part(const std::string &name) const {
3064  GMM_ASSERT1(!complex_version, "This model is a complex one");
3065  context_check();
3066  VAR_SET::iterator it = variables.find(name);
3067  GMM_ASSERT1(it != variables.end(), "Undefined variable " << name);
3068  GMM_ASSERT1(it->second.is_affine_dependent,
3069  "Only for affine dependent variables");
3070  if (act_size_to_be_done && it->second.is_fem_dofs) {
3071  if (it->second.filter != VDESCRFILTER_NO)
3072  actualize_sizes();
3073  else
3074  it->second.set_size();
3075  }
3076  for (auto &v_num : it->second.v_num_data) v_num = act_counter();
3077  return it->second.affine_real_value;
3078  }
3079 
3080  model_complex_plain_vector &
3081  model::set_complex_constant_part(const std::string &name) const {
3082  GMM_ASSERT1(complex_version, "This model is a real one");
3083  context_check();
3084  VAR_SET::iterator it = variables.find(name);
3085  GMM_ASSERT1(it!=variables.end(), "Undefined variable " << name);
3086  if (act_size_to_be_done && it->second.is_fem_dofs) {
3087  if (it->second.filter != VDESCRFILTER_NO)
3088  actualize_sizes();
3089  else
3090  it->second.set_size();
3091  }
3092  for (auto &v_num : it->second.v_num_data) v_num = act_counter();
3093  return it->second.affine_complex_value;
3094  }
3095 
3097 
3098  const brick_description &brick = bricks[ind_brick];
3099  update_brick(ind_brick, model::BUILD_ALL);
3100 
3101  brick.pbr->check_stiffness_matrix_and_rhs(*this, ind_brick, brick.tlist,
3102  brick.vlist, brick.dlist, brick.mims, brick.rmatlist,
3103  brick.rveclist[0], brick.rveclist_sym[0], brick.region);
3104  }
3105 
3106  void model::clear() {
3107  variables.clear();
3108  active_bricks.clear();
3109  valid_bricks.clear();
3110  real_dof_constraints.clear();
3111  complex_dof_constraints.clear();
3112  bricks.resize(0);
3113  rTM = model_real_sparse_matrix();
3114  cTM = model_complex_sparse_matrix();
3115  rrhs = model_real_plain_vector();
3116  crhs = model_complex_plain_vector();
3117  }
3118 
3119 
3120 
3121  void virtual_brick::full_asm_real_tangent_terms_(const model &md, size_type ind_brick,
3122  const model::varnamelist &term_list,
3123  const model::varnamelist &var_list,
3124  const model::mimlist &mim_list,
3125  model::real_matlist &rmatlist,
3126  model::real_veclist &rveclist,
3127  model::real_veclist &rveclist_sym,
3128  size_type region, build_version version) const
3129  {
3130  real_pre_assembly_in_serial(md, ind_brick, term_list, var_list, mim_list, rmatlist,
3131  rveclist, rveclist_sym, region, version);
3132  asm_real_tangent_terms(md, ind_brick, term_list, var_list, mim_list, rmatlist,
3133  rveclist, rveclist_sym, region, version);
3134  real_post_assembly_in_serial(md, ind_brick, term_list, var_list, mim_list, rmatlist,
3135  rveclist, rveclist_sym, region, version);
3136  }
3137 
3139  (const model &md, size_type s,
3140  const model::termlist& tlist,
3141  const model::varnamelist &vl,
3142  const model::varnamelist &dl,
3143  const model::mimlist &mims,
3144  model::real_matlist &matl,
3145  model::real_veclist &rvc1,
3146  model::real_veclist &rvc2,
3147  size_type rg,
3148  const scalar_type TINY) const
3149  {
3150  std::cout<<"******Verifying stiffnesses of *******"<<std::endl;
3151  std::cout<<"*** "<<brick_name()<<std::endl;
3152 
3153  //Build the index for the corresponding RHS
3154  std::map<std::string,size_type> rhs_index;
3155  for(size_type iterm=0;iterm<matl.size();iterm++)
3156  if (tlist[iterm].var1==tlist[iterm].var2) rhs_index[tlist[iterm].var1]=iterm;
3157 
3158  if (rhs_index.size()==0){
3159  GMM_WARNING0("*** cannot verify stiffness for this brick***");
3160  return;
3161  }
3162  full_asm_real_tangent_terms_(md, s, vl, dl, mims, matl, rvc1, rvc2,
3163  rg, model::BUILD_MATRIX);
3164  for(size_type iterm=0;iterm<matl.size();iterm++)
3165  {
3166 
3167  std::cout<<std::endl;
3168  std::cout<<" Stiffness["<<tlist[iterm].var1
3169  <<","<<tlist[iterm].var2<<"]:"<<std::endl;
3170  if (md.real_variable(tlist[iterm].var1).size()==0)
3171  {
3172  std::cout<<" "<<tlist[iterm].var1<<" has zero size. Skipping this term"<<std::endl;
3173  continue;
3174  }
3175  if (md.real_variable(tlist[iterm].var2).size()==0)
3176  {
3177  std::cout<<" "<<tlist[iterm].var2<<" has zero size. Skipping this term"<<std::endl;
3178  continue;
3179  }
3180 
3181  model_real_sparse_matrix SM(matl[iterm]);
3182  gmm::fill(rvc1[rhs_index[tlist[iterm].var1]], 0.0);
3183  full_asm_real_tangent_terms_(md, s, vl, dl, mims, matl, rvc1, rvc2,
3184  rg, model::BUILD_RHS);
3185  if (gmm::mat_euclidean_norm(matl[iterm])<1e-12){
3186  std::cout<<" The assembled matrix is nearly zero, skipping."<<std::endl;
3187  continue;
3188  }
3189  model_real_plain_vector RHS0(rvc1[rhs_index[tlist[iterm].var1]]);
3190 
3191  //finite difference stiffness
3192  model_real_sparse_matrix fdSM(matl[iterm].nrows(), matl[iterm].ncols());
3193  model_real_plain_vector&U = md.set_real_variable(tlist[iterm].var2);
3194  model_real_plain_vector& RHS1 = rvc1[rhs_index[tlist[iterm].var1]];
3195 
3196  scalar_type relative_tiny = gmm::vect_norminf(RHS1)*TINY;
3197  if (relative_tiny < TINY) relative_tiny = TINY;
3198 
3199  for (size_type j = 0; j < matl[iterm].ncols(); j++){
3200  U[j] += relative_tiny;
3201  gmm::fill(RHS1, 0.0);
3202  full_asm_real_tangent_terms_(md, s, vl, dl, mims, matl, rvc1, rvc2,
3203  rg, model::BUILD_RHS);
3204  for (size_type i = 0; i<matl[iterm].nrows(); i++)
3205  fdSM(i, j) = (RHS0[i] - RHS1[i]) / relative_tiny;
3206  U[j] -= relative_tiny;
3207  }
3208 
3209  model_real_sparse_matrix diffSM(matl[iterm].nrows(),matl[iterm].ncols());
3210  gmm::add(SM,gmm::scaled(fdSM,-1.0),diffSM);
3211  scalar_type norm_error_euc
3212  = gmm::mat_euclidean_norm(diffSM)/gmm::mat_euclidean_norm(SM)*100;
3213  scalar_type norm_error_1
3214  = gmm::mat_norm1(diffSM)/gmm::mat_norm1(SM)*100;
3215  scalar_type norm_error_max
3216  = gmm::mat_maxnorm(diffSM)/gmm::mat_maxnorm(SM)*100;
3217 
3218  //checking symmetry of diagonal terms
3219  scalar_type nsym_norm_error_euc=0.0;
3220  scalar_type nsym_norm_error_1=0.0;
3221  scalar_type nsym_norm_error_max=0.0;
3222  if (tlist[iterm].var1==tlist[iterm].var2){
3223  model_real_sparse_matrix diffSMtransposed(matl[iterm].nrows(),matl[iterm].ncols());
3224  gmm::add(gmm::transposed(fdSM),gmm::scaled(fdSM,-1.0),diffSMtransposed);
3225  nsym_norm_error_euc
3226  = gmm::mat_euclidean_norm(diffSMtransposed)/gmm::mat_euclidean_norm(fdSM)*100;
3227  nsym_norm_error_1
3228  = gmm::mat_norm1(diffSMtransposed)/gmm::mat_norm1(fdSM)*100;
3229  nsym_norm_error_max
3230  = gmm::mat_maxnorm(diffSMtransposed)/gmm::mat_maxnorm(fdSM)*100;
3231  }
3232 
3233  //print matrix if the size is small
3234  if(rvc1[0].size()<8){
3235  std::cout << "RHS Stiffness Matrix: \n";
3236  std::cout << "------------------------\n";
3237  for(size_type i=0; i < rvc1[iterm].size(); ++i){
3238  std::cout << "[";
3239  for(size_type j=0; j < rvc1[iterm].size(); ++j){
3240  std::cout << fdSM(i,j) << " ";
3241  }
3242  std::cout << "]\n";
3243  }
3244  std::cout << "Analytical Stiffness Matrix: \n";
3245  std::cout << "------------------------\n";
3246  for(size_type i=0; i < rvc1[iterm].size(); ++i){
3247  std::cout << "[";
3248  for(size_type j=0; j < rvc1[iterm].size(); ++j){
3249  std::cout << matl[iterm](i,j) << " ";
3250  }
3251  std::cout << "]\n";
3252  }
3253  std::cout << "Vector U: \n";
3254  std::cout << "------------------------\n";
3255  for(size_type i=0; i < rvc1[iterm].size(); ++i){
3256  std::cout << "[";
3257  std::cout << md.real_variable(tlist[iterm].var2)[i] << " ";
3258  std::cout << "]\n";
3259  }
3260  }
3261  std::cout
3262  << "\n\nfinite diff test error_norm_eucl: " << norm_error_euc << "%\n"
3263  << "finite diff test error_norm1: " << norm_error_1 << "%\n"
3264  << "finite diff test error_max_norm: " << norm_error_max << "%\n\n\n";
3265 
3266  if (tlist[iterm].var1==tlist[iterm].var2){
3267  std::cout
3268  << "Nonsymmetrical test error_norm_eucl: "<< nsym_norm_error_euc<< "%\n"
3269  << "Nonsymmetrical test error_norm1: " << nsym_norm_error_1 << "%\n"
3270  << "Nonsymmetrical test error_max_norm: " << nsym_norm_error_max << "%"
3271  << std::endl;
3272  }
3273  }
3274  }
3275 
3276  // ----------------------------------------------------------------------
3277  //
3278  //
3279  // Standard bricks
3280  //
3281  //
3282  // ----------------------------------------------------------------------
3283 
3284  // ----------------------------------------------------------------------
3285  //
3286  // Generic assembly source term brick
3287  //
3288  // ----------------------------------------------------------------------
3289 
3290  struct gen_source_term_assembly_brick : public virtual_brick {
3291 
3292  std::string expr, directvarname, directdataname;
3293  model::varnamelist vl_test1;
3294  std::string secondary_domain;
3295 
3296  void asm_real_tangent_terms(const model &md, size_type /* ib */,
3297  const model::varnamelist &,
3298  const model::varnamelist &,
3299  const model::mimlist &mims,
3300  model::real_matlist &,
3301  model::real_veclist &vecl,
3302  model::real_veclist &,
3303  size_type region,
3304  build_version) const override {
3305  GMM_ASSERT1(vecl.size() == vl_test1.size()
3306  + ((directdataname.size() == 0) ? 0 : 1), "Wrong number "
3307  "of terms for Generic source term assembly brick ");
3308  GMM_ASSERT1(mims.size() == 1, "Generic source term assembly brick "
3309  "needs one and only one mesh_im");
3310  GMM_TRACE2("Generic source term assembly");
3311 
3312  gmm::clear(vecl[0]);
3313 
3314  if (expr.size()) {
3315  const mesh_im &mim = *mims[0];
3316  mesh_region rg(region);
3317  mim.linked_mesh().intersect_with_mpi_region(rg);
3318 
3319  // reenables disabled variables
3320  ga_workspace workspace(md, ga_workspace::inherit::ALL);
3321  GMM_TRACE2(name << ": generic source term assembly");
3322  workspace.add_expression(expr, mim, rg, 1, secondary_domain);
3323  workspace.assembly(1);
3324  const auto &V=workspace.assembled_vector();
3325  for (size_type i = 0; i < vl_test1.size(); ++i) {
3326  const auto &I=workspace.interval_of_variable(vl_test1[i]);
3327  gmm::copy(gmm::sub_vector(V, I), vecl[i]);
3328  }
3329  }
3330 
3331  if (directvarname.size()) {
3332  gmm::copy(md.real_variable(directdataname), vecl.back());
3333  }
3334  }
3335 
3336  void real_post_assembly_in_serial(const model &md, size_type ib,
3337  const model::varnamelist &/* vl */,
3338  const model::varnamelist &/* dl */,
3339  const model::mimlist &/* mims */,
3340  model::real_matlist &/*matl*/,
3341  model::real_veclist &vecl,
3342  model::real_veclist &,
3343  size_type /*region*/,
3344  build_version) const override {
3345  scalar_type el = scalar_type(0);
3346  for (size_type i=0; i < vecl.size(); ++i) el += gmm::vect_norm1(vecl[i]);
3347  md.add_external_load(ib, el);
3348  }
3349 
3350  virtual std::string declare_volume_assembly_string
3351  (const model &, size_type, const model::varnamelist &,
3352  const model::varnamelist &) const {
3353  return std::string();
3354  }
3355 
3356  gen_source_term_assembly_brick(const std::string &expr_,
3357  std::string brickname,
3358  const model::varnamelist &vl_test1_,
3359  const std::string &directvarname_,
3360  const std::string &directdataname_,
3361  const std::string &secdom)
3362  : vl_test1(vl_test1_), secondary_domain(secdom) {
3363  if (brickname.size() == 0)
3364  brickname = "Generic source term assembly brick";
3365  expr = expr_;
3366  set_flags(brickname, true /* is linear*/,
3367  true /* is symmetric */, true /* is coercive */,
3368  true /* is real */, false /* is complex */,
3369  false /* compute each time */);
3370  directvarname = directvarname_; directdataname = directdataname_;
3371  }
3372 
3373  };
3374 
3375  static bool check_compatibility_vl_test(model &md,
3376  const model::varnamelist vl_test) {
3377  model::varnamelist org;
3378  for (size_type i = 0; i < vl_test.size(); ++i) {
3379  if (md.is_affine_dependent_variable(vl_test[i]))
3380  org.push_back(md.org_variable(vl_test[i]));
3381  }
3382  for (size_type i = 0; i < vl_test.size(); ++i)
3383  for (size_type j = 0; j < org.size(); ++j)
3384  if (vl_test[i].compare(org[j]) == 0) return false;
3385  return true;
3386  }
3387 
3388  size_type add_source_term_
3389  (model &md, const mesh_im &mim, const std::string &expr, size_type region,
3390  const std::string &brickname, std::string directvarname,
3391  const std::string &directdataname, bool return_if_nonlin,
3392  const std::string &secondary_domain) {
3393 
3394  ga_workspace workspace(md);
3395  size_type order = workspace.add_expression(expr, mim, region, 1,
3396  secondary_domain);
3397  GMM_ASSERT1(order <= 1, "Wrong order for a source term");
3398  model::varnamelist vl, vl_test1, vl_test2, dl;
3399  bool is_lin = workspace.used_variables(vl, vl_test1, vl_test2, dl, 1);
3400  if (!is_lin && return_if_nonlin) return size_type(-1);
3401  GMM_ASSERT1(is_lin, "Nonlinear term");
3402  GMM_ASSERT1(check_compatibility_vl_test(md, vl_test1),
3403  "This brick do not support the assembly on both an affine "
3404  "dependent variable and its original variable. "
3405  "Split the brick.");
3406 
3407  if (directdataname.size()) {
3408  vl.push_back(directvarname);
3409  dl.push_back(directdataname);
3410  } else directvarname = "";
3411 
3412  pbrick pbr = std::make_shared<gen_source_term_assembly_brick>
3413  (expr, brickname, vl_test1, directvarname, directdataname,
3414  secondary_domain);
3415  model::termlist tl;
3416 
3417  for (size_type i = 0; i < vl_test1.size(); ++i)
3418  tl.push_back(model::term_description(vl_test1[i]));
3419  if (directdataname.size())
3420  tl.push_back(model::term_description(directvarname));
3421 
3422  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
3423  }
3424 
3426  (model &md, const mesh_im &mim, const std::string &expr, size_type region,
3427  const std::string &brickname, const std::string &directvarname,
3428  const std::string &directdataname, bool return_if_nonlin) {
3429  return add_source_term_(md, mim, expr, region, brickname, directvarname,
3430  directdataname, return_if_nonlin, "");
3431  }
3432 
3434  (model &md, const mesh_im &mim, const std::string &expr, size_type region,
3435  const std::string &secondary_domain,
3436  const std::string &brickname, const std::string &directvarname,
3437  const std::string &directdataname, bool return_if_nonlin) {
3438  return add_source_term_(md, mim, expr, region, brickname, directvarname,
3439  directdataname, return_if_nonlin, secondary_domain);
3440  }
3441 
3442  // ----------------------------------------------------------------------
3443  //
3444  // Linear generic assembly brick
3445  //
3446  // ----------------------------------------------------------------------
3447 
3448  struct gen_linear_assembly_brick : public virtual_brick {
3449 
3450  std::string expr;
3451  bool is_lower_dim;
3452  model::varnamelist vl_test1, vl_test2;
3453  std::string secondary_domain;
3454 
3455  virtual void asm_real_tangent_terms(const model &md, size_type ib,
3456  const model::varnamelist &/* vl */,
3457  const model::varnamelist &dl,
3458  const model::mimlist &mims,
3459  model::real_matlist &matl,
3460  model::real_veclist &/* vecl */,
3461  model::real_veclist &,
3462  size_type region,
3463  build_version version) const {
3464  GMM_ASSERT1(matl.size() == vl_test1.size(),
3465  "Wrong number of terms for Generic linear assembly brick");
3466  GMM_ASSERT1(mims.size() == 1,
3467  "Generic linear assembly brick needs one and only one "
3468  "mesh_im");
3469  bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0);
3470  for (size_type i = 0; i < dl.size(); ++i)
3471  recompute_matrix = recompute_matrix ||
3472  md.is_var_newer_than_brick(dl[i], ib);
3473 
3474  if (recompute_matrix) {
3475  // reenables disabled variables
3476  ga_workspace workspace(md, ga_workspace::inherit::ALL);
3477  workspace.add_expression(expr, *(mims[0]), region, 2, secondary_domain);
3478  GMM_TRACE2(name << ": generic matrix assembly");
3479  workspace.assembly(2);
3480  const auto &R=workspace.assembled_matrix();
3481  for (size_type i = 0; i < vl_test1.size(); ++i) {
3482  scalar_type alpha = scalar_type(1)
3483  / ( workspace.factor_of_variable(vl_test1[i]) *
3484  workspace.factor_of_variable(vl_test2[i]));
3485  const auto &I1=workspace.interval_of_variable(vl_test1[i]),
3486  &I2=workspace.interval_of_variable(vl_test2[i]);
3487  gmm::copy(gmm::scaled(gmm::sub_matrix(R, I1, I2), alpha),
3488  matl[i]);
3489  }
3490  }
3491  }
3492 
3493  virtual std::string declare_volume_assembly_string
3494  (const model &, size_type, const model::varnamelist &,
3495  const model::varnamelist &) const {
3496  return is_lower_dim ? std::string() : expr;
3497  }
3498 
3499  gen_linear_assembly_brick(const std::string &expr_, const mesh_im &mim,
3500  bool is_sym,
3501  bool is_coer, std::string brickname,
3502  const model::varnamelist &vl_test1_,
3503  const model::varnamelist &vl_test2_,
3504  const std::string &secdom)
3505  : vl_test1(vl_test1_), vl_test2(vl_test2_), secondary_domain(secdom) {
3506  if (brickname.size() == 0) brickname = "Generic linear assembly brick";
3507  expr = expr_;
3508  is_lower_dim = mim.is_lower_dimensional();
3509  set_flags(brickname, true /* is linear*/,
3510  is_sym /* is symmetric */, is_coer /* is coercive */,
3511  true /* is real */, false /* is complex */);
3512  }
3513 
3514  };
3515 
3516  static bool check_compatibility_vl_test(model &md,
3517  const model::varnamelist vl_test1,
3518  const model::varnamelist vl_test2) {
3519  for (size_type i = 0; i < vl_test1.size(); ++i)
3520  for (size_type j = 0; j < vl_test1.size(); ++j) {
3521  bool is1 = md.is_affine_dependent_variable(vl_test1[i]);
3522  bool is2 = md.is_affine_dependent_variable(vl_test2[i]);
3523  if (is1 || is2) {
3524  const std::string &org1
3525  = is1 ? md.org_variable(vl_test1[i]) : vl_test1[i];
3526  const std::string &org2
3527  = is2 ? md.org_variable(vl_test2[i]) : vl_test2[i];
3528  bool is1_bis = md.is_affine_dependent_variable(vl_test1[j]);
3529  bool is2_bis = md.is_affine_dependent_variable(vl_test2[j]);
3530  const std::string &org1_bis = is1_bis ? md.org_variable(vl_test1[j])
3531  : vl_test1[j];
3532  const std::string &org2_bis = is2_bis ? md.org_variable(vl_test2[j])
3533  : vl_test2[j];
3534  if (org1.compare(org1_bis) == 0 && org2.compare(org2_bis))
3535  return false;
3536  }
3537  }
3538  return true;
3539  }
3540 
3541 
3542 
3543  size_type add_linear_term_
3544  (model &md, const mesh_im &mim, const std::string &expr, size_type region,
3545  bool is_sym, bool is_coercive, const std::string &brickname,
3546  bool return_if_nonlin, const std::string &secondary_domain) {
3547  // reenables disabled variables
3548  ga_workspace workspace(md, ga_workspace::inherit::ALL);
3549  size_type order = workspace.add_expression(expr, mim, region,
3550  2, secondary_domain);
3551  model::varnamelist vl, vl_test1, vl_test2, dl;
3552  bool is_lin = workspace.used_variables(vl, vl_test1, vl_test2, dl, 2);
3553 
3554  if (!is_lin && return_if_nonlin) return size_type(-1);
3555  GMM_ASSERT1(is_lin, "Nonlinear term");
3556  if (order == 0) { is_coercive = is_sym = true; }
3557 
3558  std::string const_expr= workspace.extract_constant_term(mim.linked_mesh());
3559  if (const_expr.size()) {
3560  add_source_term_
3561  (md, mim, const_expr, region, brickname+" (source term)",
3562  "", "", false, secondary_domain);
3563  }
3564 
3565  // GMM_ASSERT1(order <= 1,
3566  // "This brick does not support a second order term");
3567  GMM_ASSERT1(check_compatibility_vl_test(md, vl_test1, vl_test2),
3568  "This brick do not support the assembly on both an affine "
3569  "dependent variable and its original variable. "
3570  "Split the brick.");
3571 
3572  if (vl_test1.size()) {
3573  pbrick pbr = std::make_shared<gen_linear_assembly_brick>
3574  (expr, mim, is_sym, is_coercive, brickname, vl_test1, vl_test2,
3575  secondary_domain);
3576  model::termlist tl;
3577  for (size_type i = 0; i < vl_test1.size(); ++i)
3578  tl.push_back(model::term_description(vl_test1[i], vl_test2[i], false));
3579 
3580  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
3581  }
3582  return size_type(-1);
3583  }
3584 
3586  (model &md, const mesh_im &mim, const std::string &expr, size_type region,
3587  bool is_sym, bool is_coercive, const std::string &brickname,
3588  bool return_if_nonlin) {
3589  return add_linear_term_(md, mim, expr, region, is_sym, is_coercive,
3590  brickname, return_if_nonlin, "");
3591  }
3592 
3594  (model &md, const mesh_im &mim, const std::string &expr, size_type region,
3595  const std::string &secondary_domain, bool is_sym, bool is_coercive,
3596  const std::string &brickname, bool return_if_nonlin) {
3597  return add_linear_term_(md, mim, expr, region, is_sym, is_coercive,
3598  brickname, return_if_nonlin, secondary_domain);
3599  }
3600 
3601 
3602  // ----------------------------------------------------------------------
3603  //
3604  // Nonlinear generic assembly brick
3605  //
3606  // ----------------------------------------------------------------------
3607 
3608  struct gen_nonlinear_assembly_brick : public virtual_brick {
3609 
3610  std::string expr;
3611  bool is_lower_dim;
3612  std::string secondary_domain;
3613 
3614  virtual void real_post_assembly_in_serial(const model &md, size_type ,
3615  const model::varnamelist &,
3616  const model::varnamelist &,
3617  const model::mimlist &mims,
3618  model::real_matlist &,
3619  model::real_veclist &,
3620  model::real_veclist &,
3621  size_type region,
3622  build_version) const {
3623  GMM_ASSERT1(mims.size() == 1,
3624  "Generic linear assembly brick needs one and only one "
3625  "mesh_im");
3626  md.add_generic_expression(expr, *(mims[0]), region, secondary_domain);
3627  }
3628 
3629  virtual std::string declare_volume_assembly_string
3630  (const model &, size_type, const model::varnamelist &,
3631  const model::varnamelist &) const {
3632  return expr;
3633  }
3634 
3635 
3636  gen_nonlinear_assembly_brick(const std::string &expr_, const mesh_im &mim,
3637  bool is_sym,
3638  bool is_coer,
3639  std::string brickname,
3640  const std::string &secdom) {
3641  if (brickname.size() == 0) brickname = "Generic linear assembly brick";
3642  expr = expr_;
3643  secondary_domain = secdom;
3644  is_lower_dim = mim.is_lower_dimensional();
3645  set_flags(brickname, false /* is linear*/,
3646  is_sym /* is symmetric */, is_coer /* is coercive */,
3647  true /* is real */, false /* is complex */);
3648  }
3649 
3650  };
3651 
3652  size_type add_nonlinear_term_
3653  (model &md, const mesh_im &mim, const std::string &expr, size_type region,
3654  bool is_sym, bool is_coercive, const std::string &brickname,
3655  const std::string &secondary_domain) {
3656 
3657  ga_workspace workspace(md);
3658  size_type order = workspace.add_expression(expr, mim, region, 2,
3659  secondary_domain);
3660  GMM_ASSERT1(order < 2, "Order two test functions (Test2) are not allowed"
3661  " in assembly string for nonlinear terms");
3662  model::varnamelist vl, vl_test1, vl_test2, ddl, dl;
3663  workspace.used_variables(vl, vl_test1, vl_test2, ddl, order);
3664  for (size_type i = 0; i < ddl.size(); ++i)
3665  if (md.is_true_data(ddl[i])) dl.push_back(ddl[i]);
3666  else vl.push_back(ddl[i]);
3667  if (order == 0) { is_coercive = is_sym = true; }
3668  pbrick pbr = std::make_shared<gen_nonlinear_assembly_brick>
3669  (expr, mim, is_sym, is_coercive, brickname, secondary_domain);
3670  model::termlist tl; // No term
3671  // tl.push_back(model::term_description(true, is_sym));
3672  // TODO to be changed.
3673  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
3674  }
3675 
3676 
3678  (model &md, const mesh_im &mim, const std::string &expr, size_type region,
3679  bool is_sym, bool is_coercive, const std::string &brickname) {
3680  return add_nonlinear_term_(md, mim, expr, region, is_sym, is_coercive,
3681  brickname, "");
3682  }
3683 
3685  (model &md, const mesh_im &mim, const std::string &expr, size_type region,
3686  const std::string &secondary_domain, bool is_sym, bool is_coercive,
3687  const std::string &brickname) {
3688  return add_nonlinear_term_(md, mim, expr, region, is_sym, is_coercive,
3689  brickname, secondary_domain);
3690  }
3691 
3692 
3693  // ----------------------------------------------------------------------
3694  //
3695  // Generic elliptic brick
3696  //
3697  // ----------------------------------------------------------------------
3698 
3699  // Kept only for the complex version
3700  struct generic_elliptic_brick : public virtual_brick {
3701 
3702  virtual void asm_real_tangent_terms(const model &md, size_type /*ib*/,
3703  const model::varnamelist &vl,
3704  const model::varnamelist &dl,
3705  const model::mimlist &mims,
3706  model::real_matlist &matl,
3707  model::real_veclist &,
3708  model::real_veclist &,
3709  size_type region,
3710  build_version) const {
3711  GMM_ASSERT1(matl.size() == 1,
3712  "Generic elliptic brick has one and only one term");
3713  GMM_ASSERT1(mims.size() == 1,
3714  "Generic elliptic brick need one and only one mesh_im");
3715  GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
3716  "Wrong number of variables for generic elliptic brick");
3717 
3718  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
3719  const mesh &m = mf_u.linked_mesh();
3720  size_type N = m.dim(), Q = mf_u.get_qdim(), s = 1;
3721  const mesh_im &mim = *mims[0];
3722  const model_real_plain_vector *A = 0;
3723  const mesh_fem *mf_a = 0;
3724  mesh_region rg(region);
3725  m.intersect_with_mpi_region(rg);
3726  if (dl.size() > 0) {
3727  A = &(md.real_variable(dl[0]));
3728  mf_a = md.pmesh_fem_of_variable(dl[0]);
3729  s = gmm::vect_size(*A);
3730  if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
3731  }
3732 
3733  gmm::clear(matl[0]);
3734  GMM_TRACE2("Generic elliptic term assembly");
3735  if (s == 1) {
3736  if (mf_a) {
3737  if (Q > 1)
3739  (matl[0], mim, mf_u, *mf_a, *A, rg);
3740  else
3742  (matl[0], mim, mf_u, *mf_a, *A, rg);
3743 
3744  } else {
3745  if (Q > 1)
3747  (matl[0], mim, mf_u, rg);
3748  else
3750  (matl[0], mim, mf_u, rg);
3751  if (A) gmm::scale(matl[0], (*A)[0]);
3752  }
3753  } else if (s == N*N) {
3754  if (mf_a) {
3755  if (Q > 1)
3757  (matl[0], mim, mf_u, *mf_a, *A, rg);
3758  else
3760  (matl[0], mim, mf_u, *mf_a, *A, rg);
3761  } else {
3762  if (Q > 1)
3764  (matl[0], mim, mf_u, *A, rg);
3765  else
3767  (matl[0], mim, mf_u, *A, rg);
3768  }
3769  } else if (s == N*N*Q*Q) {
3770  if (mf_a)
3772  (matl[0], mim, mf_u, *mf_a, *A, rg);
3773  else
3775  (matl[0], mim, mf_u, *A, rg);
3776  } else
3777  GMM_ASSERT1(false, "Bad format generic elliptic brick coefficient");
3778 
3779  }
3780 
3781  virtual void real_post_assembly_in_serial(const model &md, size_type,
3782  const model::varnamelist &,
3783  const model::varnamelist &dl,
3784  const model::mimlist &/* mims */,
3785  model::real_matlist &/*matl*/,
3786  model::real_veclist &,
3787  model::real_veclist &,
3788  size_type /*region*/,
3789  build_version) const
3790  {
3791  const model_real_plain_vector *A = 0;
3792  const mesh_fem *mf_a = 0;
3793  if (dl.size() > 0)
3794  {
3795  A = &(md.real_variable(dl[0]));
3796  mf_a = md.pmesh_fem_of_variable(dl[0]);
3797  }
3798  }
3799 
3800  virtual void complex_post_assembly_in_serial(const model &md, size_type,
3801  const model::varnamelist &,
3802  const model::varnamelist &dl,
3803  const model::mimlist &/*mims*/,
3804  model::complex_matlist &/*matl*/,
3805  model::complex_veclist &,
3806  model::complex_veclist &,
3807  size_type /* region */,
3808  build_version) const
3809  {
3810  const model_real_plain_vector *A = 0;
3811  const mesh_fem *mf_a = 0;
3812  if (dl.size() > 0)
3813  {
3814  A = &(md.real_variable(dl[0]));
3815  mf_a = md.pmesh_fem_of_variable(dl[0]);
3816  }
3817  }
3818 
3819  virtual scalar_type asm_complex_tangent_terms(const model &md, size_type,
3820  const model::varnamelist &vl,
3821  const model::varnamelist &,
3822  const model::mimlist &,
3823  model::complex_matlist &matl,
3824  model::complex_veclist &,
3825  model::complex_veclist &,
3826  size_type) const {
3827  const model_complex_plain_vector &U = md.complex_variable(vl[0]);
3828  return gmm::abs(gmm::vect_hp(matl[0], U, U)) / scalar_type(2);
3829  }
3830 
3831 
3832  virtual void asm_complex_tangent_terms(const model &md, size_type,
3833  const model::varnamelist &vl,
3834  const model::varnamelist &dl,
3835  const model::mimlist &mims,
3836  model::complex_matlist &matl,
3837  model::complex_veclist &,
3838  model::complex_veclist &,
3839  size_type region,
3840  build_version) const {
3841  GMM_ASSERT1(matl.size() == 1,
3842  "Generic elliptic brick has one and only one term");
3843  GMM_ASSERT1(mims.size() == 1,
3844  "Generic elliptic brick need one and only one mesh_im");
3845  GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
3846  "Wrong number of variables for generic elliptic brick");
3847 
3848  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
3849  const mesh &m = mf_u.linked_mesh();
3850  size_type N = m.dim(), Q = mf_u.get_qdim(), s = 1;
3851  const mesh_im &mim = *mims[0];
3852  const model_real_plain_vector *A = 0;
3853  const mesh_fem *mf_a = 0;
3854  mesh_region rg(region);
3855  m.intersect_with_mpi_region(rg);
3856 
3857 
3858  if (dl.size() > 0) {
3859  A = &(md.real_variable(dl[0]));
3860  mf_a = md.pmesh_fem_of_variable(dl[0]);
3861  s = gmm::vect_size(*A);
3862  if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
3863  }
3864 
3865  gmm::clear(matl[0]);
3866  GMM_TRACE2("Generic elliptic term assembly");
3867  if (s == 1) {
3868  if (mf_a) {
3869  if (Q > 1)
3871  (matl[0], mim, mf_u, *mf_a, *A, rg);
3872  else
3874  (matl[0], mim, mf_u, *mf_a, *A, rg);
3875 
3876  } else {
3877  if (Q > 1)
3879  (gmm::real_part(matl[0]), mim, mf_u, rg);
3880  else
3882  (gmm::real_part(matl[0]), mim, mf_u, rg);
3883  if (A) gmm::scale(matl[0], (*A)[0]);
3884  }
3885  } else if (s == N*N) {
3886  if (mf_a) {
3887  if (Q > 1)
3889  (matl[0], mim, mf_u, *mf_a, *A, rg);
3890  else
3892  (matl[0], mim, mf_u, *mf_a, *A, rg);
3893  } else {
3894  if (Q > 1)
3896  (matl[0], mim, mf_u, *A, rg);
3897  else
3899  (matl[0], mim, mf_u, *A, rg);
3900  }
3901  } else if (s == N*N*Q*Q) {
3902  if (mf_a)
3904  (matl[0], mim, mf_u, *mf_a, *A, rg);
3905  else
3907  (matl[0], mim, mf_u, *A, rg);
3908  } else
3909  GMM_ASSERT1(false,
3910  "Bad format generic elliptic brick coefficient");
3911  }
3912 
3913  generic_elliptic_brick() {
3914  set_flags("Generic elliptic", true /* is linear*/,
3915  true /* is symmetric */, true /* is coercive */,
3916  true /* is real */, true /* is complex */);
3917  }
3918 
3919  };
3920 
3922  const std::string &varname,
3923  size_type region) {
3924  if (md.is_complex()) {
3925  pbrick pbr = std::make_shared<generic_elliptic_brick>();
3926  model::termlist tl;
3927  tl.push_back(model::term_description(varname, varname, true));
3928  return md.add_brick(pbr, model::varnamelist(1, varname),
3929  model::varnamelist(), tl, model::mimlist(1, &mim),
3930  region);
3931  } else {
3932  std::string test_varname
3933  = "Test_" + sup_previous_and_dot_to_varname(varname);
3934  const mesh_fem &mf_u = md.mesh_fem_of_variable(varname);
3935  size_type qdim = mf_u.get_qdim();
3936  std::string expr;
3937  if (qdim == 1)
3938  expr = "Grad_"+varname+".Grad_"+test_varname;
3939  else
3940  expr = "Grad_"+varname+":Grad_"+test_varname;
3941  return add_linear_term(md, mim, expr, region, true, true,
3942  "Laplacian", false);
3943  }
3944  }
3945 
3947  const std::string &varname,
3948  const std::string &dataname,
3949  size_type region) {
3950  if (md.is_complex()) {
3951  pbrick pbr = std::make_shared<generic_elliptic_brick>();
3952  model::termlist tl;
3953  tl.push_back(model::term_description(varname, varname, true));
3954  return md.add_brick(pbr, model::varnamelist(1, varname),
3955  model::varnamelist(1, dataname), tl,
3956  model::mimlist(1, &mim), region);
3957  } else {
3958  std::string test_varname
3959  = "Test_" + sup_previous_and_dot_to_varname(varname);
3960  const mesh_fem &mf_u = md.mesh_fem_of_variable(varname);
3961  size_type dim = mf_u.linked_mesh().dim(), qdim = mf_u.get_qdim(), qdim_data = 1;
3962  std::string expr;
3963 
3964  if (md.variable_exists(dataname)) {
3965  const mesh_fem *mf = md.pmesh_fem_of_variable(dataname);
3966  size_type n = gmm::vect_size(md.real_variable(dataname));
3967  if (mf) qdim_data = mf->get_qdim() * (n / mf->nb_dof());
3968  else qdim_data = n;
3969  }
3970 
3971  if (qdim == 1) {
3972  if (qdim_data != 1) {
3973  GMM_ASSERT1(qdim_data == gmm::sqr(dim),
3974  "Wrong data size for generic elliptic brick");
3975  expr = "((Reshape("+dataname+",meshdim,meshdim))*Grad_"+varname+").Grad_"
3976  + test_varname;
3977  } else {
3978  expr = "(("+dataname+")*Grad_"+varname+").Grad_"+test_varname;
3979  }
3980  } else {
3981  if (qdim_data != 1) {
3982  if (qdim_data == gmm::sqr(dim))
3983  expr = "((Reshape("+dataname+",meshdim,meshdim))*Grad_"+varname+"):Grad_"
3984  +test_varname;
3985  else if (qdim_data == gmm::sqr(gmm::sqr(dim)))
3986  expr = "((Reshape("+dataname+",meshdim,meshdim,meshdim,meshdim))*Grad_"
3987  +varname+"):Grad_"+test_varname;
3988  else GMM_ASSERT1(false, "Wrong data size for generic elliptic brick");
3989  } else {
3990  expr = "(("+dataname+")*Grad_"+varname+"):Grad_"+test_varname;
3991  }
3992  }
3994  (md, mim, expr, region, true, true, "Generic elliptic", true);
3995  if (ib == size_type(-1))
3996  ib = add_nonlinear_term(md, mim, expr, region, false, false,
3997  "Generic elliptic (nonlinear)");
3998  return ib;
3999  }
4000  }
4001 
4002  // ----------------------------------------------------------------------
4003  //
4004  // Source term brick
4005  //
4006  // ----------------------------------------------------------------------
4007 
4008  // Kept only for the complex version
4009  struct source_term_brick : public virtual_brick {
4010 
4011  void asm_real_tangent_terms(const model &md, size_type /*ib*/,
4012  const model::varnamelist &vl,
4013  const model::varnamelist &dl,
4014  const model::mimlist &mims,
4015  model::real_matlist &,
4016  model::real_veclist &vecl,
4017  model::real_veclist &,
4018  size_type region,
4019  build_version) const override {
4020  GMM_ASSERT1(vecl.size() == 1,
4021  "Source term brick has one and only one term");
4022  GMM_ASSERT1(mims.size() == 1,
4023  "Source term brick need one and only one mesh_im");
4024  GMM_ASSERT1(vl.size() == 1 && dl.size() > 0 && dl.size() <= 2,
4025  "Wrong number of variables for source term brick");
4026 
4027  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4028  const mesh_im &mim = *mims[0];
4029  const model_real_plain_vector &A = md.real_variable(dl[0]);
4030  const mesh_fem *mf_data = md.pmesh_fem_of_variable(dl[0]);
4031  mesh_region rg(region);
4032  mim.linked_mesh().intersect_with_mpi_region(rg);
4033 
4034  size_type s = gmm::vect_size(A);
4035  if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4036 
4037  GMM_ASSERT1(mf_u.get_qdim() == s,
4038  dl[0] << ": bad format of source term data. "
4039  "Detected dimension is " << s << " should be "
4040  << size_type(mf_u.get_qdim()));
4041 
4042  GMM_TRACE2("Source term assembly");
4043  if (mf_data)
4044  asm_source_term(vecl[0], mim, mf_u, *mf_data, A, rg);
4045  else
4046  asm_homogeneous_source_term(vecl[0], mim, mf_u, A, rg);
4047 
4048  if (dl.size() > 1) gmm::add(md.real_variable(dl[1]), vecl[0]);
4049  }
4050 
4051  void real_post_assembly_in_serial(const model &md, size_type ib,
4052  const model::varnamelist &/* vl */,
4053  const model::varnamelist &/* dl */,
4054  const model::mimlist &/* mims */,
4055  model::real_matlist &/*matl*/,
4056  model::real_veclist &vecl,
4057  model::real_veclist &,
4058  size_type /*region*/,
4059  build_version) const override
4060  {
4061  md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4062  }
4063 
4064 
4065  void asm_complex_tangent_terms(const model &md, size_type /*ib*/,
4066  const model::varnamelist &vl,
4067  const model::varnamelist &dl,
4068  const model::mimlist &mims,
4069  model::complex_matlist &,
4070  model::complex_veclist &vecl,
4071  model::complex_veclist &,
4072  size_type region,
4073  build_version) const override {
4074  GMM_ASSERT1(vecl.size() == 1,
4075  "Source term brick has one and only one term");
4076  GMM_ASSERT1(mims.size() == 1,
4077  "Source term brick need one and only one mesh_im");
4078  GMM_ASSERT1(vl.size() == 1 && dl.size() > 0 && dl.size() <= 2,
4079  "Wrong number of variables for source term brick");
4080 
4081  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4082  const mesh_im &mim = *mims[0];
4083  const model_complex_plain_vector &A = md.complex_variable(dl[0]);
4084  const mesh_fem *mf_data = md.pmesh_fem_of_variable(dl[0]);
4085  mesh_region rg(region);
4086  mim.linked_mesh().intersect_with_mpi_region(rg);
4087 
4088  size_type s = gmm::vect_size(A);
4089  if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4090 
4091  GMM_ASSERT1(mf_u.get_qdim() == s, "Bad format of source term data");
4092 
4093  GMM_TRACE2("Source term assembly");
4094  if (mf_data)
4095  asm_source_term(vecl[0], mim, mf_u, *mf_data, A, rg);
4096  else
4097  asm_homogeneous_source_term(vecl[0], mim, mf_u, A, rg);
4098 
4099  if (dl.size() > 1) gmm::add(md.complex_variable(dl[1]), vecl[0]);
4100  }
4101 
4102  void complex_post_assembly_in_serial(const model &md,
4103  size_type ib,
4104  const model::varnamelist &,
4105  const model::varnamelist &,
4106  const model::mimlist &,
4107  model::complex_matlist &,
4108  model::complex_veclist &vecl,
4109  model::complex_veclist &,
4110  size_type, build_version) const override
4111  {
4112  md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4113  }
4114 
4115 
4116 
4117  source_term_brick() {
4118  set_flags("Source term", true /* is linear*/,
4119  true /* is symmetric */, true /* is coercive */,
4120  true /* is real */, true /* is complex */,
4121  false /* compute each time */);
4122  }
4123 
4124 
4125  };
4126 
4128  const std::string &varname,
4129  const std::string &dataexpr,
4130  size_type region,
4131  const std::string &directdataname) {
4132  if (md.is_complex()) {
4133  pbrick pbr = std::make_shared<source_term_brick>();
4134  model::termlist tl;
4135  tl.push_back(model::term_description(varname));
4136  model::varnamelist vdata(1, dataexpr);
4137  if (directdataname.size()) vdata.push_back(directdataname);
4138  return md.add_brick(pbr, model::varnamelist(1, varname),
4139  vdata, tl, model::mimlist(1, &mim), region);
4140  } else {
4141  std::string test_varname
4142  = "Test_" + sup_previous_and_dot_to_varname(varname);
4143  const mesh_fem &mf_u = md.mesh_fem_of_variable(varname);
4144  size_type qdim = mf_u.get_qdim();
4145  std::string expr;
4146  if (qdim == 1)
4147  expr = "("+dataexpr+")*"+test_varname;
4148  else
4149  expr = "("+dataexpr+")."+test_varname;
4150  size_type ib = add_source_term_generic_assembly_brick
4151  (md, mim, expr, region, "Source term", varname, directdataname, true);
4152  if (ib == size_type(-1)) {
4153  ib = add_nonlinear_term(md, mim, "-("+expr+")", region, false, false,
4154  "Source term (nonlinear)");
4155  if (directdataname.size())
4156  add_source_term_generic_assembly_brick
4157  (md, mim, "", region, "Source term", varname, directdataname);
4158  }
4159  return ib;
4160  }
4161  }
4162 
4163  // ----------------------------------------------------------------------
4164  //
4165  // Normal source term brick
4166  //
4167  // ----------------------------------------------------------------------
4168 
4169  struct normal_source_term_brick : public virtual_brick {
4170 
4171  void asm_real_tangent_terms(const model &md, size_type /* ib */,
4172  const model::varnamelist &vl,
4173  const model::varnamelist &dl,
4174  const model::mimlist &mims,
4175  model::real_matlist &,
4176  model::real_veclist &vecl,
4177  model::real_veclist &,
4178  size_type region,
4179  build_version) const override {
4180  GMM_ASSERT1(vecl.size() == 1,
4181  "Source term brick has one and only one term");
4182  GMM_ASSERT1(mims.size() == 1,
4183  "Source term brick need one and only one mesh_im");
4184  GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
4185  "Wrong number of variables for source term brick");
4186 
4187  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4188  const mesh_im &mim = *mims[0];
4189  const model_real_plain_vector &A = md.real_variable(dl[0]);
4190  const mesh_fem *mf_data = md.pmesh_fem_of_variable(dl[0]);
4191  mesh_region rg(region);
4192  mim.linked_mesh().intersect_with_mpi_region(rg);
4193 
4194  size_type s = gmm::vect_size(A), N = mf_u.linked_mesh().dim();
4195  if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4196 
4197  GMM_ASSERT1(mf_u.get_qdim()*N == s,
4198  dl[0] << ": bad format of normal source term data. "
4199  "Detected dimension is " << s << " should be "
4200  << size_type(mf_u.get_qdim()*N));
4201 
4202  GMM_TRACE2("source term assembly");
4203  if (mf_data)
4204  asm_normal_source_term(vecl[0], mim, mf_u, *mf_data, A, rg);
4205  else
4206  asm_homogeneous_normal_source_term(vecl[0], mim, mf_u, A, rg);
4207  }
4208 
4209  void real_post_assembly_in_serial(const model &md, size_type ib,
4210  const model::varnamelist &/* vl */,
4211  const model::varnamelist &/* dl */,
4212  const model::mimlist &/* mims */,
4213  model::real_matlist &/*matl*/,
4214  model::real_veclist &vecl,
4215  model::real_veclist &,
4216  size_type /*region*/,
4217  build_version) const override {
4218  md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4219  }
4220 
4221 
4222  virtual void asm_complex_tangent_terms(const model &md, size_type /* ib */,
4223  const model::varnamelist &vl,
4224  const model::varnamelist &dl,
4225  const model::mimlist &mims,
4226  model::complex_matlist &,
4227  model::complex_veclist &vecl,
4228  model::complex_veclist &,
4229  size_type region,
4230  build_version) const {
4231  GMM_ASSERT1(vecl.size() == 1,
4232  "Source term brick has one and only one term");
4233  GMM_ASSERT1(mims.size() == 1,
4234  "Source term brick need one and only one mesh_im");
4235  GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
4236  "Wrong number of variables for source term brick");
4237 
4238  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4239  const mesh_im &mim = *mims[0];
4240  const model_complex_plain_vector &A = md.complex_variable(dl[0]);
4241  const mesh_fem *mf_data = md.pmesh_fem_of_variable(dl[0]);
4242  mesh_region rg(region);
4243  mim.linked_mesh().intersect_with_mpi_region(rg);
4244 
4245  size_type s = gmm::vect_size(A), N = mf_u.linked_mesh().dim();
4246  if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4247 
4248  GMM_ASSERT1(s == mf_u.get_qdim()*N, "Bad format of source term data");
4249 
4250  GMM_TRACE2("source term assembly");
4251  if (mf_data)
4252  asm_normal_source_term(vecl[0], mim, mf_u, *mf_data, A, rg);
4253  else
4254  asm_homogeneous_normal_source_term(vecl[0], mim, mf_u, A, rg);
4255  }
4256 
4257  void complex_post_assembly_in_serial(const model &md, size_type ib,
4258  const model::varnamelist &/* vl */,
4259  const model::varnamelist &/* dl */,
4260  const model::mimlist &/* mims */,
4261  model::complex_matlist &/*matl*/,
4262  model::complex_veclist &vecl,
4263  model::complex_veclist &,
4264  size_type /*region*/,
4265  build_version) const override {
4266  md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4267  }
4268 
4269  normal_source_term_brick() {
4270  set_flags("Normal source term", true /* is linear*/,
4271  true /* is symmetric */, true /* is coercive */,
4272  true /* is real */, true /* is complex */,
4273  false /* compute each time */);
4274  }
4275 
4276 
4277  };
4278 
4280  const std::string &varname,
4281  const std::string &dataexpr,
4282  size_type region) {
4283  if (md.is_complex()) {
4284  pbrick pbr = std::make_shared<normal_source_term_brick>();
4285  model::termlist tl;
4286  tl.push_back(model::term_description(varname));
4287  model::varnamelist vdata(1, dataexpr);
4288  return md.add_brick(pbr, model::varnamelist(1, varname),
4289  vdata, tl, model::mimlist(1, &mim), region);
4290  } else {
4291  std::string test_varname
4292  = "Test_" + sup_previous_and_dot_to_varname(varname);
4293  const mesh_fem &mf_u = md.mesh_fem_of_variable(varname);
4294  size_type qdim = mf_u.get_qdim();
4295  std::string expr;
4296  if (qdim == 1)
4297  expr = "(("+dataexpr+").Normal)*"+test_varname;
4298  else
4299  expr = "(Reshape("+dataexpr+",qdim("+varname
4300  + "),meshdim)*Normal)."+test_varname;
4301  return add_source_term_generic_assembly_brick
4302  (md, mim, expr, region, "Source term");
4303  }
4304  }
4305 
4306 
4307  // ----------------------------------------------------------------------
4308  //
4309  // Dirichlet condition brick
4310  //
4311  // ----------------------------------------------------------------------
4312  // Two variables : with multipliers
4313  // One variable : penalization
4314 
4315  struct Dirichlet_condition_brick : public virtual_brick {
4316 
4317  bool H_version; // The version hu = r for vector fields.
4318  bool normal_component; // Dirichlet on normal component for vector field.
4319  const mesh_fem *mf_mult_;
4324 
4325  virtual void asm_real_tangent_terms(const model &md, size_type ib,
4326  const model::varnamelist &vl,
4327  const model::varnamelist &dl,
4328  const model::mimlist &mims,
4329  model::real_matlist &matl,
4330  model::real_veclist &vecl,
4331  model::real_veclist &,
4332  size_type region,
4333  build_version version) const {
4334  GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
4335  "Dirichlet condition brick has one and only one term");
4336  GMM_ASSERT1(mims.size() == 1,
4337  "Dirichlet condition brick need one and only one mesh_im");
4338  GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() <= 3,
4339  "Wrong number of variables for Dirichlet condition brick");
4340 
4341  model_real_sparse_matrix& rB = rB_th;
4342  model_real_plain_vector& rV = rV_th;
4343 
4344  bool penalized = (vl.size() == 1);
4345  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4346  const mesh_fem &mf_mult = penalized ? (mf_mult_ ? *mf_mult_ : mf_u)
4347  : md.mesh_fem_of_variable(vl[1]);
4348  const mesh_im &mim = *mims[0];
4349  const model_real_plain_vector *A = 0, *COEFF = 0, *H = 0;
4350  const mesh_fem *mf_data = 0, *mf_H = 0;
4351  bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
4352  || (penalized && md.is_var_newer_than_brick(dl[0], ib));
4353 
4354  if (penalized) {
4355  COEFF = &(md.real_variable(dl[0]));
4356  GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
4357  "Data for coefficient should be a scalar");
4358  }
4359 
4360  size_type s = 0, ind = (penalized ? 1 : 0);
4361  if (dl.size() > ind) {
4362  A = &(md.real_variable(dl[ind]));
4363  mf_data = md.pmesh_fem_of_variable(dl[ind]);
4364  s = gmm::vect_size(*A);
4365  if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4366  size_type ss = ((normal_component) ? 1 : mf_u.get_qdim());
4367  GMM_ASSERT1(s == ss, dl[ind] << ": bad format of "
4368  "Dirichlet data. Detected dimension is " << s
4369  << " should be " << ss);
4370  }
4371 
4372  if (dl.size() > ind + 1) {
4373  GMM_ASSERT1(H_version,
4374  "Wrong number of data for Dirichlet condition brick");
4375  H = &(md.real_variable(dl[ind+1]));
4376  mf_H = md.pmesh_fem_of_variable(dl[ind+1]);
4377  s = gmm::vect_size(*A);
4378  if (mf_H) {
4379  s = s * mf_H->get_qdim() / mf_H->nb_dof();
4380  GMM_ASSERT1(mf_H->get_qdim() == 1, "Implemented only for mf_H "
4381  "a scalar finite element method");
4382  }
4383  GMM_ASSERT1(s = gmm::sqr(mf_u.get_qdim()),
4384  dl[ind+1] << ": bad format of Dirichlet data. "
4385  "Detected dimension is " << s << " should be "
4386  << size_type(gmm::sqr(mf_u.get_qdim())));
4387  }
4388 
4389  mesh_region rg(region);
4390  mim.linked_mesh().intersect_with_mpi_region(rg);
4391 
4392  if (recompute_matrix) {
4393  model_real_sparse_matrix *B = &(matl[0]);
4394  if (penalized && (&mf_mult != &mf_u)) {
4395  gmm::resize(rB, mf_mult.nb_dof(), mf_u.nb_dof());
4396  gmm::clear(rB);
4397  B = &rB;
4398  } else {
4399  gmm::clear(matl[0]);
4400  }
4401  GMM_TRACE2("Mass term assembly for Dirichlet condition");
4402  if (H_version || normal_component) {
4403  ga_workspace workspace;
4404  gmm::sub_interval Imult(0, mf_mult.nb_dof()), Iu(0, mf_u.nb_dof());
4405  base_vector u(mf_u.nb_dof());
4406  base_vector mult(mf_mult.nb_dof());
4407  workspace.add_fem_variable("u", mf_u, Iu, u);
4408  workspace.add_fem_variable("mult", mf_mult, Imult, mult);
4409  auto expression = std::string{};
4410  if (H_version){
4411  if (mf_H) workspace.add_fem_constant("A", *mf_H, *H);
4412  else workspace.add_fixed_size_constant("A", *H);
4413  expression = (mf_u.get_qdim() == 1) ?
4414  "Test_mult . (A . Test2_u)"
4415  :
4416  "Test_mult. (Reshape(A, qdim(u), qdim(u)) . Test2_u)";
4417  } else if (normal_component) {
4418  expression = "Test_mult . (Test2_u . Normal)";
4419  }
4420  workspace.add_expression(expression, mim, rg);
4421  workspace.set_assembled_matrix(*B);
4422  workspace.assembly(2);
4423  } else {
4424  asm_mass_matrix(*B, mim, mf_mult, mf_u, rg);
4425  }
4426 
4427  if (penalized && (&mf_mult != &mf_u)) {
4428  GMM_ASSERT1(MPI_IS_MASTER(), "Sorry, the penalized option "
4429  "filtered by a multiplier space is not parallelized");
4430  gmm::mult(gmm::transposed(rB), rB, matl[0]);
4431  gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
4432  } else if (penalized) {
4433  gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
4434  }
4435  }
4436 
4437  if (dl.size() > ind) {
4438  GMM_TRACE2("Source term assembly for Dirichlet condition");
4439 
4440  if (penalized && (&mf_mult != &mf_u)) {
4441  gmm::resize(rV, mf_mult.nb_dof());
4442  gmm::clear(rV);
4443  if (mf_data)
4444  asm_source_term(rV, mim, mf_mult, *mf_data, *A, rg);
4445  else
4446  asm_homogeneous_source_term(rV, mim, mf_mult, *A, rg);
4447  } else {
4448  if (mf_data)
4449  asm_source_term(vecl[0], mim, mf_mult, *mf_data, *A, rg);
4450  else
4451  asm_homogeneous_source_term(vecl[0], mim, mf_mult, *A, rg);
4452  }
4453 
4454  if (penalized && (&mf_mult != &mf_u)) {
4455  gmm::mult(gmm::transposed(rB), rV, vecl[0]);
4456  gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
4457  rV = model_real_plain_vector();
4458  } else if (penalized)
4459  gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
4460  }
4461 
4462  }
4463 
4464  void real_post_assembly_in_serial(const model &md, size_type ib,
4465  const model::varnamelist &/* vl */,
4466  const model::varnamelist &/* dl */,
4467  const model::mimlist &/* mims */,
4468  model::real_matlist &/*matl*/,
4469  model::real_veclist &vecl,
4470  model::real_veclist &,
4471  size_type /*region*/,
4472  build_version) const override {
4473  md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4474  }
4475 
4476  virtual void asm_complex_tangent_terms(const model &md, size_type ib,
4477  const model::varnamelist &vl,
4478  const model::varnamelist &dl,
4479  const model::mimlist &mims,
4480  model::complex_matlist &matl,
4481  model::complex_veclist &vecl,
4482  model::complex_veclist &,
4483  size_type region,
4484  build_version version) const {
4485  GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
4486  "Dirichlet condition brick has one and only one term");
4487  GMM_ASSERT1(mims.size() == 1,
4488  "Dirichlet condition brick need one and only one mesh_im");
4489  GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() <= 3,
4490  "Wrong number of variables for Dirichlet condition brick");
4491 
4492  model_complex_sparse_matrix& cB = cB_th;
4493  model_complex_plain_vector& cV = cV_th;
4494 
4495  bool penalized = (vl.size() == 1);
4496  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4497  const mesh_fem &mf_mult = penalized ? (mf_mult_ ? *mf_mult_ : mf_u)
4498  : md.mesh_fem_of_variable(vl[1]);
4499  const mesh_im &mim = *mims[0];
4500  const model_complex_plain_vector *A = 0, *COEFF = 0, *H = 0;
4501  const mesh_fem *mf_data = 0, *mf_H = 0;
4502  bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
4503  || (penalized && md.is_var_newer_than_brick(dl[0], ib));
4504 
4505  if (penalized) {
4506  COEFF = &(md.complex_variable(dl[0]));
4507  GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
4508  "Data for coefficient should be a scalar");
4509  }
4510 
4511  size_type s = 0, ind = (penalized ? 1 : 0);
4512  if (dl.size() > ind) {
4513  A = &(md.complex_variable(dl[ind]));
4514  mf_data = md.pmesh_fem_of_variable(dl[ind]);
4515  s = gmm::vect_size(*A);
4516  if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
4517  size_type ss = s * ((normal_component) ? mf_u.linked_mesh().dim() : 1);
4518  GMM_ASSERT1(mf_u.get_qdim() == ss,
4519  dl[ind] << ": bad format of Dirichlet data. "
4520  "Detected dimension is " << ss << " should be "
4521  << size_type(mf_u.get_qdim()));
4522  }
4523 
4524  if (dl.size() > ind + 1) {
4525  GMM_ASSERT1(H_version,
4526  "Wrong number of data for Dirichlet condition brick");
4527  H = &(md.complex_variable(dl[ind+1]));
4528  mf_H = md.pmesh_fem_of_variable(dl[ind+1]);
4529  s = gmm::vect_size(*A);
4530  if (mf_H) {
4531  s = s * mf_H->get_qdim() / mf_H->nb_dof();
4532  GMM_ASSERT1(mf_H->get_qdim() == 1, "Implemented only for mf_H "
4533  "a scalar finite element method");
4534  }
4535  GMM_ASSERT1(s = gmm::sqr(mf_u.get_qdim()),
4536  dl[ind+1] << ": bad format of Dirichlet data. "
4537  "Detected dimension is " << s << " should be "
4538  << size_type(gmm::sqr(mf_u.get_qdim())));
4539  }
4540 
4541  mesh_region rg(region);
4542  mim.linked_mesh().intersect_with_mpi_region(rg);
4543 
4544  if (recompute_matrix) {
4545  model_complex_sparse_matrix *B = &(matl[0]);
4546  if (penalized && (&mf_mult != &mf_u)) {
4547  gmm::resize(cB, mf_mult.nb_dof(), mf_u.nb_dof());
4548  gmm::clear(cB);
4549  B = &cB;
4550  } else {
4551  gmm::clear(matl[0]);
4552  }
4553  GMM_TRACE2("Mass term assembly for Dirichlet condition");
4554  if (H_version) {
4555  if (mf_u.get_qdim() == 1)
4556  asm_real_or_complex_1_param_mat(*B, mim, mf_mult, mf_H, *H, rg,
4557  "(A*Test_u).Test2_u");
4558  else
4559  asm_real_or_complex_1_param_mat(*B, mim, mf_mult, mf_H, *H, rg,
4560  "(Reshape(A,qdim(u),qdim(u))*Test2_u).Test_u");
4561  // if (mf_H)
4562  // asm_real_or_complex_1_param
4563  // (*B, mim, mf_mult, *mf_H, *H, rg, (mf_u.get_qdim() == 1) ?
4564  // "F=data(#2);"
4565  // "M(#1,#3)+=sym(comp(Base(#1).Base(#3).Base(#2))(:,:,i).F(i))"
4566  // : "F=data(qdim(#1),qdim(#1),#2);"
4567  // "M(#1,#3)+=sym(comp(vBase(#1).vBase(#3).Base(#2))(:,i,:,j,k).F(i,j,k));", &mf_u);
4568  // else
4569  // asm_real_or_complex_1_param
4570  // (*B, mim, mf_mult, mf_u, *H, rg, (mf_u.get_qdim() == 1) ?
4571  // "F=data(1);"
4572  // "M(#1,#2)+=sym(comp(Base(#1).Base(#2)).F(1))"
4573  // : "F=data(qdim(#1),qdim(#1));"
4574  // "M(#1,#2)+=sym(comp(vBase(#1).vBase(#2))(:,i,:,j).F(i,j));");
4575  }
4576  else if (normal_component) {
4577  ga_workspace workspace;
4578  gmm::sub_interval Imult(0, mf_mult.nb_dof()), Iu(0, mf_u.nb_dof());
4579  base_vector mult(mf_mult.nb_dof()), u(mf_u.nb_dof());
4580  workspace.add_fem_variable("mult", mf_mult, Imult, mult);
4581  workspace.add_fem_variable("u", mf_u, Iu, u);
4582  workspace.add_expression("Test_mult.(Test2_u.Normal)", mim, rg);
4583  model_real_sparse_matrix BB(mf_mult.nb_dof(), mf_u.nb_dof());
4584  workspace.set_assembled_matrix(BB);
4585  workspace.assembly(2);
4586  gmm::add(BB, *B);
4587 
4588  // generic_assembly assem;
4589  // if (mf_mult.get_qdim() == 1)
4590  // assem.set("M(#2,#1)+=comp(Base(#2).vBase(#1).Normal())(:,:,i,i);");
4591  // else
4592  // assem.set("M(#2,#1)+=comp(vBase(#2).mBase(#1).Normal())(:,i,:,i,j,j);");
4593  // assem.push_mi(mim);
4594  // assem.push_mf(mf_u);
4595  // assem.push_mf(mf_mult);
4596  // assem.push_mat(gmm::real_part(*B));
4597  // assem.assembly(rg);
4598  } else {
4599  asm_mass_matrix(*B, mim, mf_mult, mf_u, rg);
4600  }
4601  if (penalized && (&mf_mult != &mf_u)) {
4602  gmm::mult(gmm::transposed(cB), cB, matl[0]);
4603  gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
4604  } else if (penalized) {
4605  gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
4606  }
4607  }
4608 
4609  if (dl.size() > ind) {
4610  GMM_TRACE2("Source term assembly for Dirichlet condition");
4611 
4612  if (penalized && (&mf_mult != &mf_u)) {
4613  gmm::resize(cV, mf_mult.nb_dof());
4614  gmm::clear(cV);
4615  if (mf_data)
4616  asm_source_term(cV, mim, mf_mult, *mf_data, *A, rg);
4617  else
4618  asm_homogeneous_source_term(cV, mim, mf_mult, *A, rg);
4619  } else {
4620  if (mf_data)
4621  asm_source_term(vecl[0], mim, mf_mult, *mf_data, *A, rg);
4622  else
4623  asm_homogeneous_source_term(vecl[0], mim, mf_mult, *A, rg);
4624  }
4625 
4626  if (penalized && (&mf_mult != &mf_u)) {
4627  gmm::mult(gmm::transposed(cB), cV, vecl[0]);
4628  gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
4629  cV = model_complex_plain_vector();
4630  } else if (penalized)
4631  gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
4632  }
4633  }
4634 
4635  void complex_post_assembly_in_serial(const model &md, size_type ib,
4636  const model::varnamelist &/* vl */,
4637  const model::varnamelist &/* dl */,
4638  const model::mimlist &/* mims */,
4639  model::complex_matlist &/*matl*/,
4640  model::complex_veclist &vecl,
4641  model::complex_veclist &,
4642  size_type /*region*/,
4643  build_version) const override {
4644  md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
4645  }
4646 
4647 
4648  virtual std::string declare_volume_assembly_string
4649  (const model &, size_type, const model::varnamelist &,
4650  const model::varnamelist &) const {
4651  return std::string();
4652  }
4653 
4654  Dirichlet_condition_brick(bool penalized, bool H_version_,
4655  bool normal_component_,
4656  const mesh_fem *mf_mult__ = 0) {
4657  mf_mult_ = mf_mult__;
4658  H_version = H_version_;
4659  normal_component = normal_component_;
4660  GMM_ASSERT1(!(H_version && normal_component), "Bad Dirichlet version");
4661  set_flags(penalized ? "Dirichlet with penalization brick"
4662  : "Dirichlet with multipliers brick",
4663  true /* is linear*/,
4664  true /* is symmetric */, penalized /* is coercive */,
4665  true /* is real */, true /* is complex */,
4666  false /* compute each time */);
4667  }
4668  };
4669 
4671  (model &md, const mesh_im &mim, const std::string &varname,
4672  const std::string &multname, size_type region,
4673  const std::string &dataname) {
4674  pbrick pbr = std::make_shared<Dirichlet_condition_brick>(false,false,false);
4675  model::termlist tl;
4676  tl.push_back(model::term_description(multname, varname, true));
4677  model::varnamelist vl(1, varname);
4678  vl.push_back(multname);
4679  model::varnamelist dl;
4680  if (dataname.size()) dl.push_back(dataname);
4681  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4682  }
4683 
4685  (model &md, const mesh_im &mim, const std::string &varname,
4686  const mesh_fem &mf_mult, size_type region,
4687  const std::string &dataname) {
4688  std::string multname = md.new_name("mult_on_" + varname);
4689  md.add_multiplier(multname, mf_mult, varname);
4691  (md, mim, varname, multname, region, dataname);
4692  }
4693 
4695  (model &md, const mesh_im &mim, const std::string &varname,
4696  dim_type degree, size_type region,
4697  const std::string &dataname) {
4698  const mesh_fem &mf_u = md.mesh_fem_of_variable(varname);
4699  const mesh_fem &mf_mult = classical_mesh_fem(mf_u.linked_mesh(),
4700  degree, mf_u.get_qdim());
4702  (md, mim, varname, mf_mult, region, dataname);
4703  }
4704 
4705  const std::string &mult_varname_Dirichlet(model &md, size_type ind_brick) {
4706  return md.varname_of_brick(ind_brick, 1);
4707  }
4708 
4710  (model &md, const mesh_im &mim, const std::string &varname,
4711  scalar_type penalisation_coeff, size_type region,
4712  const std::string &dataname, const mesh_fem *mf_mult) {
4713  std::string coeffname = md.new_name("penalization_on_" + varname);
4714  md.add_fixed_size_data(coeffname, 1);
4715  if (md.is_complex())
4716  md.set_complex_variable(coeffname)[0] = penalisation_coeff;
4717  else
4718  md.set_real_variable(coeffname)[0] = penalisation_coeff;
4719  pbrick pbr = std::make_shared<Dirichlet_condition_brick>
4720  (true, false, false, mf_mult);
4721  model::termlist tl;
4722  tl.push_back(model::term_description(varname, varname, true));
4723  model::varnamelist vl(1, varname);
4724  model::varnamelist dl(1, coeffname);
4725  if (dataname.size()) dl.push_back(dataname);
4726  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4727  }
4728 
4730  (model &md, const mesh_im &mim, const std::string &varname,
4731  const std::string &multname, size_type region,
4732  const std::string &dataname) {
4733  pbrick pbr = std::make_shared<Dirichlet_condition_brick>(false,false,true);
4734  model::termlist tl;
4735  tl.push_back(model::term_description(multname, varname, true));
4736  model::varnamelist vl(1, varname);
4737  vl.push_back(multname);
4738  model::varnamelist dl;
4739  if (dataname.size()) dl.push_back(dataname);
4740  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4741  }
4742 
4744  (model &md, const mesh_im &mim, const std::string &varname,
4745  const mesh_fem &mf_mult, size_type region,
4746  const std::string &dataname) {
4747  std::string multname = md.new_name("mult_on_" + varname);
4748  md.add_multiplier(multname, mf_mult, varname);
4750  (md, mim, varname, multname, region, dataname);
4751  }
4752 
4754  (model &md, const mesh_im &mim, const std::string &varname,
4755  dim_type degree, size_type region,
4756  const std::string &dataname) {
4757  const mesh_fem &mf_u = md.mesh_fem_of_variable(varname);
4758  const mesh_fem &mf_mult = classical_mesh_fem(mf_u.linked_mesh(),degree, 1);
4760  (md, mim, varname, mf_mult, region, dataname);
4761  }
4762 
4764  (model &md, const mesh_im &mim, const std::string &varname,
4765  scalar_type penalisation_coeff, size_type region,
4766  const std::string &dataname, const mesh_fem *mf_mult) {
4767  std::string coeffname = md.new_name("penalization_on_" + varname);
4768  md.add_fixed_size_data(coeffname, 1);
4769  if (md.is_complex())
4770  md.set_complex_variable(coeffname)[0] = penalisation_coeff;
4771  else
4772  md.set_real_variable(coeffname)[0] = penalisation_coeff;
4773  pbrick pbr = std::make_shared<Dirichlet_condition_brick>
4774  (true, false, true, mf_mult);
4775  model::termlist tl;
4776  tl.push_back(model::term_description(varname, varname, true));
4777  model::varnamelist vl(1, varname);
4778  model::varnamelist dl(1, coeffname);
4779  if (dataname.size()) dl.push_back(dataname);
4780  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4781  }
4782 
4783 
4785  (model &md, const mesh_im &mim, const std::string &varname,
4786  const std::string &multname, size_type region,
4787  const std::string &dataname, const std::string &Hname) {
4788  pbrick pbr = std::make_shared<Dirichlet_condition_brick>(false,true,false);
4789  model::termlist tl;
4790  tl.push_back(model::term_description(multname, varname, true));
4791  model::varnamelist vl(1, varname);
4792  vl.push_back(multname);
4793  model::varnamelist dl;
4794  dl.push_back(dataname);
4795  dl.push_back(Hname);
4796  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4797  }
4798 
4800  (model &md, const mesh_im &mim, const std::string &varname,
4801  const mesh_fem &mf_mult, size_type region,
4802  const std::string &dataname, const std::string &Hname) {
4803  std::string multname = md.new_name("mult_on_" + varname);
4804  md.add_multiplier(multname, mf_mult, varname);
4806  (md, mim, varname, multname, region, dataname, Hname);
4807  }
4808 
4810  (model &md, const mesh_im &mim, const std::string &varname,
4811  dim_type degree, size_type region,
4812  const std::string &dataname, const std::string &Hname) {
4813  const mesh_fem &mf_u = md.mesh_fem_of_variable(varname);
4814  const mesh_fem &mf_mult = classical_mesh_fem(mf_u.linked_mesh(),
4815  degree, mf_u.get_qdim());
4817  (md, mim, varname, mf_mult, region, dataname, Hname);
4818  }
4819 
4821  (model &md, const mesh_im &mim, const std::string &varname,
4822  scalar_type penalisation_coeff, size_type region,
4823  const std::string &dataname, const std::string &Hname,
4824  const mesh_fem *mf_mult) {
4825  std::string coeffname = md.new_name("penalization_on_" + varname);
4826  md.add_fixed_size_data(coeffname, 1);
4827  if (md.is_complex())
4828  md.set_complex_variable(coeffname)[0] = penalisation_coeff;
4829  else
4830  md.set_real_variable(coeffname)[0] = penalisation_coeff;
4831  pbrick pbr = std::make_shared<Dirichlet_condition_brick>
4832  (true, true, false, mf_mult);
4833  model::termlist tl;
4834  tl.push_back(model::term_description(varname, varname, true));
4835  model::varnamelist vl(1, varname);
4836  model::varnamelist dl(1, coeffname);
4837  dl.push_back(dataname);
4838  dl.push_back(Hname);
4839  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
4840  }
4841 
4843  scalar_type penalisation_coeff) {
4844  const std::string &coeffname = md.dataname_of_brick(ind_brick, 0);
4845  if (!md.is_complex()) {
4846  model_real_plain_vector &d = md.set_real_variable(coeffname);
4847  GMM_ASSERT1(gmm::vect_size(d)==1,
4848  "Wrong coefficient size, may be not a Dirichlet brick "
4849  "with penalization");
4850  d[0] = penalisation_coeff;
4851  }
4852  else {
4853  model_complex_plain_vector &d = md.set_complex_variable(coeffname);
4854  GMM_ASSERT1(gmm::vect_size(d)==1,
4855  "Wrong coefficient size, may be not a Dirichlet brick "
4856  "with penalization");
4857  d[0] = penalisation_coeff;
4858  }
4859  }
4860 
4861  // ----------------------------------------------------------------------
4862  //
4863  // Dirichlet condition brick with simplification
4864  //
4865  // ----------------------------------------------------------------------
4866 
4867  struct simplification_Dirichlet_condition_brick : public virtual_brick {
4868 
4869  virtual void asm_real_tangent_terms(const model &md, size_type /*ib*/,
4870  const model::varnamelist &vl,
4871  const model::varnamelist &dl,
4872  const model::mimlist &mims,
4873  model::real_matlist &matl,
4874  model::real_veclist &vecl,
4875  model::real_veclist &,
4876  size_type region,
4877  build_version /*version*/) const {
4878  if (MPI_IS_MASTER()) {
4879 
4880  GMM_ASSERT1(vecl.size() == 0 && matl.size() == 0,
4881  "Dirichlet condition brick by simplification has no term");
4882  GMM_ASSERT1(mims.size() == 0,
4883  "Dirichlet condition brick by simplification need no "
4884  "mesh_im");
4885  GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
4886  "Wrong number of variables for Dirichlet condition brick "
4887  "by simplification");
4888 
4889  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4890  const model_real_plain_vector *A = 0;
4891  const mesh_fem *mf_data = 0;
4892  size_type s = 0;
4893 
4894  if (dl.size() == 1) {
4895  A = &(md.real_variable(dl[0]));
4896  mf_data = md.pmesh_fem_of_variable(dl[0]);
4897 
4898  if (mf_data) {
4899  GMM_ASSERT1(mf_data == &mf_u, "Sorry, for this brick, the data has"
4900  " to be defined on the same f.e.m. as the unknown");
4901  } else {
4902  s = gmm::vect_size(*A);
4903  GMM_ASSERT1(mf_u.get_qdim() == s, ": bad format of "
4904  "Dirichlet data. Detected dimension is " << s
4905  << " should be " << size_type(mf_u.get_qdim()));
4906  }
4907  }
4908 
4909  mesh_region rg(region);
4910  // mf_u.linked_mesh().intersect_with_mpi_region(rg); // Not distributed
4911 
4912  if (mf_u.get_qdim() > 1 || (!mf_data && A)) {
4913  for (mr_visitor i(rg, mf_u.linked_mesh()); !i.finished(); ++i) {
4914  pfem pf = mf_u.fem_of_element(i.cv());
4915  if (pf) {
4916  GMM_ASSERT1(pf->target_dim() == 1,
4917  "Intrinsically vectorial fems are not allowed");
4918  GMM_ASSERT1(mf_data || pf->is_lagrange(), "Constant Dirichlet "
4919  "data allowed for lagrange fems only");
4920  }
4921  }
4922  }
4923 
4924  dal::bit_vector dofs = mf_u.dof_on_region(rg);
4925 
4926  if (A && !mf_data) {
4927  GMM_ASSERT1(dofs.card() % s == 0, "Problem with dof vectorization");
4928  }
4929 
4930  for (dal::bv_visitor i(dofs); !i.finished(); ++i) {
4931  scalar_type val(0);
4932  if (A) val = (mf_data ? (*A)[i] : (*A)[i%s]);
4933  md.add_real_dof_constraint(vl[0], i, val);
4934  }
4935  }
4936  }
4937 
4938  virtual void asm_complex_tangent_terms(const model &md, size_type /*ib*/,
4939  const model::varnamelist &vl,
4940  const model::varnamelist &dl,
4941  const model::mimlist &mims,
4942  model::complex_matlist &matl,
4943  model::complex_veclist &vecl,
4944  model::complex_veclist &,
4945  size_type region,
4946  build_version /*version*/) const {
4947  if (MPI_IS_MASTER()) {
4948  GMM_ASSERT1(vecl.size() == 0 && matl.size() == 0,
4949  "Dirichlet condition brick by simplification has no term");
4950  GMM_ASSERT1(mims.size() == 0,
4951  "Dirichlet condition brick by simplification need no "
4952  "mesh_im");
4953  GMM_ASSERT1(vl.size() == 1 && dl.size() <= 1,
4954  "Wrong number of variables for Dirichlet condition brick "
4955  "by simplification");
4956 
4957  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
4958  const model_complex_plain_vector *A = 0;
4959  const mesh_fem *mf_data = 0;
4960  size_type s = 0;
4961 
4962  if (dl.size() == 1) {
4963  A = &(md.complex_variable(dl[0]));
4964  mf_data = md.pmesh_fem_of_variable(dl[0]);
4965 
4966  if (mf_data) {
4967  GMM_ASSERT1(mf_data == &mf_u, "Sorry, for this brick, the data has"
4968  " to be define on the same f.e.m. than the unknown");
4969  } else {
4970  s = gmm::vect_size(*A);
4971  GMM_ASSERT1(mf_u.get_qdim() == s, ": bad format of "
4972  "Dirichlet data. Detected dimension is " << s
4973  << " should be " << size_type(mf_u.get_qdim()));
4974  }
4975  }
4976 
4977  mesh_region rg(region);
4978  // mf_u.linked_mesh().intersect_with_mpi_region(rg); // Not distributed
4979 
4980  if (mf_u.get_qdim() > 1 || (!mf_data && A)) {
4981  for (mr_visitor i(rg, mf_u.linked_mesh()); !i.finished(); ++i) {
4982  pfem pf = mf_u.fem_of_element(i.cv());
4983  if (pf) {
4984  GMM_ASSERT1(pf->target_dim() == 1,
4985  "Intrinsically vectorial fems are not allowed");
4986  GMM_ASSERT1(mf_data || pf->is_lagrange(), "Constant Dirichlet "
4987  "data allowed for lagrange fems only");
4988  }
4989  }
4990  }
4991 
4992  dal::bit_vector dofs = mf_u.dof_on_region(rg);
4993 
4994  if (A && !mf_data) {
4995  GMM_ASSERT1(dofs.card() % s == 0, "Problem with dof vectorization");
4996  }
4997 
4998  for (dal::bv_visitor i(dofs); !i.finished(); ++i) {
4999  complex_type val(0);
5000  if (A) val = (mf_data ? (*A)[i] : (*A)[i%s]);
5001  md.add_complex_dof_constraint(vl[0], i, val);
5002  }
5003  }
5004  }
5005 
5006 
5007  virtual std::string declare_volume_assembly_string
5008  (const model &, size_type, const model::varnamelist &,
5009  const model::varnamelist &) const {
5010  return std::string();
5011  }
5012 
5013  simplification_Dirichlet_condition_brick() {
5014  set_flags("Dirichlet with simplification brick",
5015  true /* is linear*/,
5016  true /* is symmetric */, true /* is coercive */,
5017  true /* is real */, true /* is complex */,
5018  true /* compute each time */);
5019  }
5020  };
5021 
5023  (model &md, const std::string &varname,
5024  size_type region, const std::string &dataname) {
5025  pbrick pbr = std::make_shared<simplification_Dirichlet_condition_brick>();
5026  model::termlist tl;
5027  model::varnamelist vl(1, varname);
5028  model::varnamelist dl;
5029  if (dataname.size()) dl.push_back(dataname);
5030  return md.add_brick(pbr, vl, dl, tl, model::mimlist(), region);
5031  }
5032 
5033  // ----------------------------------------------------------------------
5034  //
5035  // Dirichlet condition brick with Nitsche's method
5036  //
5037  // ----------------------------------------------------------------------
5038 
5040  (model &md, const mesh_im &mim, const std::string &varname,
5041  const std::string &Neumannterm,
5042  const std::string &datagamma0, size_type region, scalar_type theta_,
5043  const std::string &datag) {
5044  std::string theta = std::to_string(theta_);
5045  // reenables disabled variables
5046  ga_workspace workspace(md, ga_workspace::inherit::ALL);
5047  size_type order = workspace.add_expression(Neumannterm, mim, region, 1);
5048  GMM_ASSERT1(order == 0, "Wrong expression of the Neumann term");
5049  bool is_lin = workspace.is_linear(1);
5050 
5051  std::string condition = "("+varname + (datag.size() ? "-("+datag+"))":")");
5052  std::string gamma = "(("+datagamma0+")*element_size)";
5053  std::string r = "(1/"+gamma+")";
5054  std::string expr = "("+r+"*"+condition+"-("+Neumannterm+")).Test_"+varname;
5055  if (theta_ != scalar_type(0)) {
5056  std::string derivative_Neumann = workspace.extract_order1_term(varname);
5057  if (derivative_Neumann.size())
5058  expr+="-"+theta+"*"+condition+".("+derivative_Neumann+")";
5059  }
5060 
5061  // cout << "Nitsche expression : " << expr << endl;
5062  // cout << "is_lin : " << int(is_lin) << endl;
5063 
5064  if (is_lin) {
5065  return add_linear_term(md, mim, expr, region, false, false,
5066  "Dirichlet condition with Nitsche's method");
5067  } else {
5068  return add_nonlinear_term(md, mim, expr, region, false, false,
5069  "Dirichlet condition with Nitsche's method");
5070  }
5071  }
5072 
5074  (model &md, const mesh_im &mim, const std::string &varname,
5075  const std::string &Neumannterm,
5076  const std::string &datagamma0, size_type region, scalar_type theta_,
5077  const std::string &datag) {
5078  std::string theta = std::to_string(theta_);
5079  // reenables disabled variables
5080  ga_workspace workspace(md, ga_workspace::inherit::ALL);
5081  size_type order = workspace.add_expression(Neumannterm, mim, region, 1);
5082  GMM_ASSERT1(order == 0, "Wrong expression of the Neumann term");
5083  bool is_lin = workspace.is_linear(1);
5084 
5085  std::string condition = "("+varname+".Normal"
5086  + (datag.size() ? "-("+datag+"))":")");
5087  std::string gamma = "(("+datagamma0+")*element_size)";
5088  std::string r = "(1/"+gamma+")";
5089  std::string expr = "("+r+"*"+condition+"-Normal.("+Neumannterm
5090  +"))*(Normal.Test_"+varname+")";
5091  if (theta_ != scalar_type(0)) {
5092  std::string derivative_Neumann = workspace.extract_order1_term(varname);
5093  if (derivative_Neumann.size())
5094  expr+="-"+theta+"*"+condition+"*Normal.("
5095  +derivative_Neumann+")";
5096  }
5097  if (is_lin) {
5098  return add_linear_term(md, mim, expr, region, false, false,
5099  "Dirichlet condition with Nitsche's method");
5100  } else {
5101  return add_nonlinear_term(md, mim, expr, region, false, false,
5102  "Dirichlet condition with Nitsche's method");
5103  }
5104  }
5105 
5107  (model &md, const mesh_im &mim, const std::string &varname,
5108  const std::string &Neumannterm,
5109  const std::string &datagamma0, size_type region, scalar_type theta_,
5110  const std::string &datag, const std::string &dataH) {
5111  std::string theta = std::to_string(theta_);
5112  // reenables disabled variables
5113  ga_workspace workspace(md, ga_workspace::inherit::ALL);
5114  size_type order = workspace.add_expression(Neumannterm, mim, region, 1);
5115  GMM_ASSERT1(order == 0, "Wrong expression of the Neumann term");
5116  bool is_lin = workspace.is_linear(1);
5117 
5118  std::string condition = "(("+dataH+")*"+varname
5119  + (datag.size() ? "-("+datag+"))":")");
5120  std::string gamma = "(("+datagamma0+")*element_size)";
5121  std::string r = "(1/"+gamma+")";
5122  std::string expr = "("+r+"*"+condition+"-("+dataH+")*("+Neumannterm
5123  +"))*(("+dataH+")*Test_"+varname+")";
5124  if (theta_ != scalar_type(0)) {
5125  std::string derivative_Neumann = workspace.extract_order1_term(varname);
5126  if (derivative_Neumann.size())
5127  expr+="-"+theta+"*"+condition+"*(("+dataH+")*("
5128  +derivative_Neumann+"))";
5129  }
5130  if (is_lin) {
5131  return add_linear_term(md, mim, expr, region, false, false,
5132  "Dirichlet condition with Nitsche's method");
5133  } else {
5134  return add_nonlinear_term(md, mim, expr, region, false, false,
5135  "Dirichlet condition with Nitsche's method");
5136  }
5137  }
5138 
5139  // ----------------------------------------------------------------------
5140  //
5141  // Pointwise constraints brick
5142  //
5143  // ----------------------------------------------------------------------
5144  // Two variables : with multipliers
5145  // One variable : penalization
5146 
5147  struct pointwise_constraints_brick : public virtual_brick {
5148 
5149  mutable gmm::row_matrix<model_real_sparse_vector> rB;
5150  mutable gmm::row_matrix<model_complex_sparse_vector> cB;
5151 
5152  virtual void real_pre_assembly_in_serial(const model &md, size_type ib,
5153  const model::varnamelist &vl,
5154  const model::varnamelist &dl,
5155  const model::mimlist &mims,
5156  model::real_matlist &matl,
5157  model::real_veclist &vecl,
5158  model::real_veclist &/*rvecl*/,
5159  size_type,
5160  build_version version) const {
5161  if (MPI_IS_MASTER()) {
5162 
5163  GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5164  "Pointwize constraints brick has only one term");
5165  GMM_ASSERT1(mims.size() == 0,
5166  "Pointwize constraints brick does not need a mesh_im");
5167  GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2,
5168  "Wrong number of variables for pointwize constraints brick");
5169  bool penalized = (vl.size() == 1);
5170  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5171  dim_type N = mf_u.linked_mesh().dim(), Q = mf_u.get_qdim(), ind_pt = 0;
5172  size_type dlsize = size_type((penalized ? 1 : 0) + 1 + (Q > 1 ? 1 : 0));
5173  GMM_ASSERT1(dl.size() == dlsize || dl.size() == dlsize+1,
5174  "Wrong number of data for pointwize constraints brick");
5175 
5176 
5177  const model_real_plain_vector *COEFF = 0;
5178  if (penalized) {
5179  COEFF = &(md.real_variable(dl[0]));
5180  ind_pt = 1;
5181  GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
5182  "Data for coefficient should be a scalar");
5183  }
5184 
5185  const model_real_plain_vector &PT = md.real_variable(dl[ind_pt]);
5186  size_type nb_co = gmm::vect_size(PT) / N;
5187 
5188  dim_type ind_unitv = dim_type((Q > 1) ? ind_pt+1 : 0);
5189  const model_real_plain_vector &unitv =md.real_variable(dl[ind_unitv]);
5190  GMM_ASSERT1((!ind_unitv || gmm::vect_size(unitv) == nb_co * Q),
5191  "Wrong size for vector of unit vectors");
5192 
5193  dim_type ind_rhs = dim_type((Q > 1) ? ind_pt+2 : ind_pt+1);
5194  if (dl.size() < size_type(ind_rhs + 1)) ind_rhs = 0;
5195  const model_real_plain_vector &rhs = md.real_variable(dl[ind_rhs]);
5196  GMM_ASSERT1((!ind_rhs || gmm::vect_size(rhs) == nb_co),
5197  "Wrong size for vector of rhs");
5198 
5199  bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
5200  || (penalized && (md.is_var_newer_than_brick(dl[ind_pt], ib)
5201  || md.is_var_newer_than_brick(dl[ind_unitv], ib)
5202  || md.is_var_newer_than_brick(dl[ind_rhs], ib)));
5203 
5204  if (recompute_matrix) {
5205  gmm::row_matrix<model_real_sparse_vector> BB(nb_co*Q, mf_u.nb_dof());
5206  gmm::clear(rB); gmm::resize(rB, nb_co, mf_u.nb_dof());
5207 
5208  dal::bit_vector dof_untouched;
5209  getfem::mesh_trans_inv mti(mf_u.linked_mesh());
5210  base_node pt(N);
5211  for (size_type i = 0; i < nb_co; ++i) {
5212  gmm::copy(gmm::sub_vector(PT, gmm::sub_interval(i*N, N)), pt);
5213  mti.add_point(pt);
5214  }
5215  gmm::row_matrix<model_real_sparse_vector> &BBB = ((Q > 1) ? BB : rB);
5216  model_real_plain_vector vv;
5217  interpolation(mf_u, mti, vv, vv, BBB, 1, 1, &dof_untouched);
5218  GMM_ASSERT1(dof_untouched.card() == 0,
5219  "Pointwize constraints : some of the points are outside "
5220  "the mesh: " << dof_untouched);
5221 
5222  if (Q > 1) {
5223  for (size_type i = 0; i < nb_co; ++i)
5224  for (size_type q = 0; q < Q; ++q)
5225  gmm::add(gmm::scaled(gmm::mat_row(BB, i*Q+q), unitv[i*Q+q]),
5226  gmm::mat_row(rB, i));
5227  }
5228  if (penalized) {
5229  gmm::mult(gmm::transposed(rB), rB, matl[0]);
5230  gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
5231  } else
5232  gmm::copy(rB, matl[0]);
5233  }
5234 
5235  if (ind_rhs) {
5236  if (penalized) {
5237  gmm::mult(gmm::transposed(rB), rhs, vecl[0]);
5238  gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
5239  }
5240  else gmm::copy(rhs, vecl[0]);
5241  }
5242  else gmm::clear(vecl[0]);
5243  }
5244  }
5245 
5246  virtual void complex_pre_assembly_in_serial(const model &md, size_type ib,
5247  const model::varnamelist &vl,
5248  const model::varnamelist &dl,
5249  const model::mimlist &mims,
5250  model::complex_matlist &matl,
5251  model::complex_veclist &vecl,
5252  model::complex_veclist &,
5253  size_type,
5254  build_version version) const {
5255  if (MPI_IS_MASTER()) {
5256  GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5257  "Pointwize constraints brick only one term");
5258  GMM_ASSERT1(mims.size() == 0,
5259  "Pointwize constraints brick does not need a mesh_im");
5260  GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2,
5261  "Wrong number of variables for pointwize constraints brick");
5262  bool penalized = (vl.size() == 1);
5263  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5264  dim_type N = mf_u.linked_mesh().dim(), Q = mf_u.get_qdim(), ind_pt = 0;
5265  size_type dlsize = size_type((penalized ? 1 : 0) + 1 + (Q > 1 ? 1 :0));
5266  GMM_ASSERT1(dl.size() == dlsize || dl.size() == dlsize+1,
5267  "Wrong number of data for pointwize constraints brick");
5268 
5269 
5270  const model_complex_plain_vector *COEFF = 0;
5271  if (penalized) {
5272  COEFF = &(md.complex_variable(dl[0]));
5273  ind_pt = 1;
5274  GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
5275  "Data for coefficient should be a scalar");
5276  }
5277 
5278  const model_complex_plain_vector &PT = md.complex_variable(dl[ind_pt]);
5279  size_type nb_co = gmm::vect_size(PT) / N;
5280 
5281  dim_type ind_unitv = dim_type((Q > 1) ? ind_pt+1 : 0);
5282  const model_complex_plain_vector &unitv
5283  = md.complex_variable(dl[ind_unitv]);
5284  GMM_ASSERT1((!ind_unitv || gmm::vect_size(unitv) == nb_co * Q),
5285  "Wrong size for vector of unit vectors");
5286 
5287  dim_type ind_rhs = dim_type((Q > 1) ? ind_pt+2 : ind_pt+1);
5288  if (dl.size() < size_type(ind_rhs + 1)) ind_rhs = 0;
5289  const model_complex_plain_vector &rhs
5290  = md.complex_variable(dl[ind_rhs]);
5291  GMM_ASSERT1((!ind_rhs || gmm::vect_size(rhs) == nb_co),
5292  "Wrong size for vector of rhs");
5293 
5294  bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
5295  || (penalized && (md.is_var_newer_than_brick(dl[ind_pt], ib)
5296  || md.is_var_newer_than_brick(dl[ind_unitv], ib)
5297  || md.is_var_newer_than_brick(dl[ind_rhs], ib)));
5298 
5299  if (recompute_matrix) {
5300  gmm::row_matrix<model_complex_sparse_vector>
5301  BB(nb_co*Q,mf_u.nb_dof());
5302  gmm::clear(cB); gmm::resize(cB, nb_co, mf_u.nb_dof());
5303  dal::bit_vector dof_untouched;
5304  getfem::mesh_trans_inv mti(mf_u.linked_mesh());
5305  base_node pt(N);
5306  for (size_type i = 0; i < nb_co; ++i) {
5307  gmm::copy(gmm::real_part
5308  (gmm::sub_vector(PT, gmm::sub_interval(i*N, N))), pt);
5309  mti.add_point(pt);
5310  }
5311  gmm::row_matrix<model_complex_sparse_vector> &BBB = ((Q > 1) ? BB :cB);
5312  model_complex_plain_vector vv;
5313  interpolation(mf_u, mti, vv, vv, BBB, 1, 1, &dof_untouched);
5314  GMM_ASSERT1(dof_untouched.card() == 0,
5315  "Pointwize constraints : some of the points are outside "
5316  "the mesh: " << dof_untouched);
5317 
5318  if (Q > 1) {
5319  for (size_type i = 0; i < nb_co; ++i)
5320  for (size_type q = 0; q < Q; ++q)
5321  gmm::add(gmm::scaled(gmm::mat_row(BB, i*Q+q), unitv[i*Q+q]),
5322  gmm::mat_row(cB, i));
5323  }
5324 
5325  if (penalized) {
5326  gmm::mult(gmm::transposed(cB), cB, matl[0]);
5327  gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
5328  } else
5329  gmm::copy(cB, matl[0]);
5330  }
5331 
5332 
5333  if (ind_rhs) {
5334  if (penalized) {
5335  gmm::mult(gmm::transposed(cB), rhs, vecl[0]);
5336  gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
5337  }
5338  else gmm::copy(rhs, vecl[0]);
5339  }
5340  else gmm::clear(vecl[0]);
5341  }
5342  }
5343 
5344  virtual std::string declare_volume_assembly_string
5345  (const model &, size_type, const model::varnamelist &,
5346  const model::varnamelist &) const {
5347  return std::string();
5348  }
5349 
5350  pointwise_constraints_brick(bool penalized) {
5351  set_flags(penalized ? "Pointwise cosntraints with penalization brick"
5352  : "Pointwise cosntraints with multipliers brick",
5353  true /* is linear*/,
5354  true /* is symmetric */, penalized /* is coercive */,
5355  true /* is real */, true /* is complex */,
5356  false /* compute each time */);
5357  }
5358  };
5359 
5360 
5362  (model &md, const std::string &varname,
5363  scalar_type penalisation_coeff, const std::string &dataname_pt,
5364  const std::string &dataname_unitv, const std::string &dataname_val) {
5365  std::string coeffname = md.new_name("penalization_on_" + varname);
5366  md.add_fixed_size_data(coeffname, 1);
5367  if (md.is_complex())
5368  md.set_complex_variable(coeffname)[0] = penalisation_coeff;
5369  else
5370  md.set_real_variable(coeffname)[0] = penalisation_coeff;
5371  pbrick pbr = std::make_shared<pointwise_constraints_brick>(true);
5372  model::termlist tl;
5373  tl.push_back(model::term_description(varname, varname, true));
5374  model::varnamelist vl(1, varname);
5375  model::varnamelist dl(1, coeffname);
5376  dl.push_back(dataname_pt);
5377  const mesh_fem &mf_u = md.mesh_fem_of_variable(varname);
5378  if (mf_u.get_qdim() > 1) dl.push_back(dataname_unitv);
5379  if (dataname_val.size() > 0) dl.push_back(dataname_val);
5380  return md.add_brick(pbr, vl, dl, tl, model::mimlist(), size_type(-1));
5381  }
5382 
5384  (model &md, const std::string &varname,
5385  const std::string &multname, const std::string &dataname_pt,
5386  const std::string &dataname_unitv, const std::string &dataname_val) {
5387  pbrick pbr = std::make_shared<pointwise_constraints_brick>(false);
5388  model::termlist tl;
5389  tl.push_back(model::term_description(multname, varname, true));
5390  model::varnamelist vl(1, varname);
5391  vl.push_back(multname);
5392  model::varnamelist dl(1, dataname_pt);
5393  const mesh_fem &mf_u = md.mesh_fem_of_variable(varname);
5394  if (mf_u.get_qdim() > 1) dl.push_back(dataname_unitv);
5395  if (dataname_val.size() > 0) dl.push_back(dataname_val);
5396  return md.add_brick(pbr, vl, dl, tl, model::mimlist(), size_type(-1));
5397  }
5398 
5400  (model &md, const std::string &varname, const std::string &dataname_pt,
5401  const std::string &dataname_unitv, const std::string &dataname_val) {
5402  std::string multname = md.new_name("mult_on_" + varname);
5403  const mesh_fem &mf_u = md.mesh_fem_of_variable(varname);
5404  size_type nb_co =
5405  ((md.is_complex()) ? gmm::vect_size(md.complex_variable(dataname_pt))
5406  : gmm::vect_size(md.real_variable(dataname_pt)))
5407  / mf_u.linked_mesh().dim();
5408  md.add_fixed_size_variable(multname, nb_co);
5410  (md, varname, multname, dataname_pt, dataname_unitv, dataname_val);
5411  }
5412 
5413 
5414  // ----------------------------------------------------------------------
5415  //
5416  // Helmholtz brick
5417  //
5418  // ----------------------------------------------------------------------
5419 
5420  struct Helmholtz_brick : public virtual_brick {
5421 
5422  virtual void asm_real_tangent_terms(const model &md, size_type,
5423  const model::varnamelist &vl,
5424  const model::varnamelist &dl,
5425  const model::mimlist &mims,
5426  model::real_matlist &matl,
5427  model::real_veclist &,
5428  model::real_veclist &,
5429  size_type region,
5430  build_version) const {
5431  GMM_ASSERT1(matl.size() == 1,
5432  "Helmholtz brick has one and only one term");
5433  GMM_ASSERT1(mims.size() == 1,
5434  "Helmholtz brick need one and only one mesh_im");
5435  GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
5436  "Wrong number of variables for Helmholtz brick");
5437 
5438  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5439  const mesh &m = mf_u.linked_mesh();
5440  size_type Q = mf_u.get_qdim(), s = 1;
5441  GMM_ASSERT1(Q == 1, "Helmholtz brick is only for scalar field, sorry.");
5442  const mesh_im &mim = *mims[0];
5443  const mesh_fem *mf_a = 0;
5444  mesh_region rg(region);
5445  m.intersect_with_mpi_region(rg);
5446  const model_real_plain_vector *A = &(md.real_variable(dl[0]));
5447  mf_a = md.pmesh_fem_of_variable(dl[0]);
5448  s = gmm::vect_size(*A);
5449  if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
5450 
5451  if (s == 1) {
5452  GMM_TRACE2("Stiffness matrix assembly for Helmholtz problem");
5453  gmm::clear(matl[0]);
5454  model_real_plain_vector A2(gmm::vect_size(*A));
5455  for (size_type i=0; i < gmm::vect_size(*A); ++i) // Not valid for
5456  A2[i] = gmm::sqr((*A)[i]); // non lagrangian fem ...
5457  if (mf_a)
5458  asm_Helmholtz(matl[0], mim, mf_u, *mf_a, A2, rg);
5459  else
5460  asm_homogeneous_Helmholtz(matl[0], mim, mf_u, A2, rg);
5461  } else
5462  GMM_ASSERT1(false, "Bad format Helmholtz brick coefficient");
5463  }
5464 
5465  virtual void asm_complex_tangent_terms(const model &md, size_type,
5466  const model::varnamelist &vl,
5467  const model::varnamelist &dl,
5468  const model::mimlist &mims,
5469  model::complex_matlist &matl,
5470  model::complex_veclist &,
5471  model::complex_veclist &,
5472  size_type region,
5473  build_version) const {
5474  GMM_ASSERT1(matl.size() == 1,
5475  "Helmholtz brick has one and only one term");
5476  GMM_ASSERT1(mims.size() == 1,
5477  "Helmholtz brick need one and only one mesh_im");
5478  GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
5479  "Wrong number of variables for Helmholtz brick");
5480 
5481  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5482  const mesh &m = mf_u.linked_mesh();
5483  size_type Q = mf_u.get_qdim(), s = 1;
5484  GMM_ASSERT1(Q == 1, "Helmholtz brick is only for scalar field, sorry.");
5485  const mesh_im &mim = *mims[0];
5486  const mesh_fem *mf_a = 0;
5487  mesh_region rg(region);
5488  m.intersect_with_mpi_region(rg);
5489  const model_complex_plain_vector *A = &(md.complex_variable(dl[0]));
5490  mf_a = md.pmesh_fem_of_variable(dl[0]);
5491  s = gmm::vect_size(*A);
5492  if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
5493 
5494  if (s == 1) {
5495  GMM_TRACE2("Stiffness matrix assembly for Helmholtz problem");
5496  gmm::clear(matl[0]);
5497  model_complex_plain_vector A2(gmm::vect_size(*A));
5498  for (size_type i=0; i < gmm::vect_size(*A); ++i) // Not valid for
5499  A2[i] = gmm::sqr((*A)[i]); // non lagrangian fem ...
5500  if (mf_a)
5501  asm_Helmholtz(matl[0], mim, mf_u, *mf_a, A2, rg);
5502  else
5503  asm_homogeneous_Helmholtz(matl[0], mim, mf_u, A2, rg);
5504  } else
5505  GMM_ASSERT1(false, "Bad format Helmholtz brick coefficient");
5506  }
5507 
5508  Helmholtz_brick() {
5509  set_flags("Helmholtz", true /* is linear*/,
5510  true /* is symmetric */, true /* is coercive */,
5511  true /* is real */, true /* is complex */);
5512  }
5513 
5514  };
5515 
5517  const std::string &varname,
5518  const std::string &dataexpr,
5519  size_type region) {
5520  if (md.is_complex()) {
5521  pbrick pbr = std::make_shared<Helmholtz_brick>();
5522  model::termlist tl;
5523  tl.push_back(model::term_description(varname, varname, true));
5524  return md.add_brick(pbr, model::varnamelist(1, varname),
5525  model::varnamelist(1, dataexpr), tl,
5526  model::mimlist(1, &mim), region);
5527  } else {
5528  std::string test_varname
5529  = "Test_" + sup_previous_and_dot_to_varname(varname);
5530  std::string expr = "Grad_"+varname+".Grad_"+test_varname
5531  +" + sqr("+dataexpr+")*"+varname+"*"+test_varname;
5532 
5533  size_type ib = add_linear_term(md, mim, expr, region, true, true,
5534  "Helmholtz", true);
5535  if (ib == size_type(-1))
5536  ib = add_nonlinear_term(md, mim, expr, region, false, false,
5537  "Helmholtz (nonlinear)");
5538  return ib;
5539  }
5540  }
5541 
5542 
5543 
5544  // ----------------------------------------------------------------------
5545  //
5546  // Fourier-Robin brick
5547  //
5548  // ----------------------------------------------------------------------
5549 
5550  struct Fourier_Robin_brick : public virtual_brick {
5551 
5552  virtual void asm_real_tangent_terms(const model &md, size_type,
5553  const model::varnamelist &vl,
5554  const model::varnamelist &dl,
5555  const model::mimlist &mims,
5556  model::real_matlist &matl,
5557  model::real_veclist &,
5558  model::real_veclist &,
5559  size_type region,
5560  build_version) const {
5561  GMM_ASSERT1(matl.size() == 1,
5562  "Fourier-Robin brick has one and only one term");
5563  GMM_ASSERT1(mims.size() == 1,
5564  "Fourier-Robin brick need one and only one mesh_im");
5565  GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
5566  "Wrong number of variables for Fourier-Robin brick");
5567 
5568  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5569  const mesh &m = mf_u.linked_mesh();
5570  size_type Q = mf_u.get_qdim(), s = 1;
5571  const mesh_im &mim = *mims[0];
5572  const mesh_fem *mf_a = 0;
5573  mesh_region rg(region);
5574  m.intersect_with_mpi_region(rg);
5575  const model_real_plain_vector *A = &(md.real_variable(dl[0]));
5576  mf_a = md.pmesh_fem_of_variable(dl[0]);
5577  s = gmm::vect_size(*A);
5578  if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
5579  GMM_ASSERT1(s == Q*Q,
5580  "Bad format Fourier-Robin brick coefficient");
5581 
5582  GMM_TRACE2("Fourier-Robin term assembly");
5583  gmm::clear(matl[0]);
5584  if (mf_a)
5585  asm_qu_term(matl[0], mim, mf_u, *mf_a, *A, rg);
5586  else
5587  asm_homogeneous_qu_term(matl[0], mim, mf_u, *A, rg);
5588  }
5589 
5590  virtual void asm_complex_tangent_terms(const model &md, size_type,
5591  const model::varnamelist &vl,
5592  const model::varnamelist &dl,
5593  const model::mimlist &mims,
5594  model::complex_matlist &matl,
5595  model::complex_veclist &,
5596  model::complex_veclist &,
5597  size_type region,
5598  build_version) const {
5599  GMM_ASSERT1(matl.size() == 1,
5600  "Fourier-Robin brick has one and only one term");
5601  GMM_ASSERT1(mims.size() == 1,
5602  "Fourier-Robin brick need one and only one mesh_im");
5603  GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
5604  "Wrong number of variables for Fourier-Robin brick");
5605 
5606  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
5607  const mesh &m = mf_u.linked_mesh();
5608  size_type Q = mf_u.get_qdim(), s = 1;
5609  const mesh_im &mim = *mims[0];
5610  const mesh_fem *mf_a = 0;
5611  mesh_region rg(region);
5612  m.intersect_with_mpi_region(rg);
5613  const model_complex_plain_vector *A = &(md.complex_variable(dl[0]));
5614  mf_a = md.pmesh_fem_of_variable(dl[0]);
5615  s = gmm::vect_size(*A);
5616  if (mf_a) s = s * mf_a->get_qdim() / mf_a->nb_dof();
5617  GMM_ASSERT1(s == Q*Q,
5618  "Bad format Fourier-Robin brick coefficient");
5619 
5620  GMM_TRACE2("Fourier-Robin term assembly");
5621  gmm::clear(matl[0]);
5622  if (mf_a)
5623  asm_qu_term(matl[0], mim, mf_u, *mf_a, *A, rg);
5624  else
5625  asm_homogeneous_qu_term(matl[0], mim, mf_u, *A, rg);
5626  }
5627 
5628  Fourier_Robin_brick() {
5629  set_flags("Fourier Robin condition", true /* is linear*/,
5630  true /* is symmetric */, true /* is coercive */,
5631  true /* is real */, true /* is complex */,
5632  false /* compute each time */);
5633  }
5634 
5635  };
5636 
5638  const std::string &varname,
5639  const std::string &dataexpr,
5640  size_type region) {
5641  if (md.is_complex()) {
5642  pbrick pbr = std::make_shared<Fourier_Robin_brick>();
5643  model::termlist tl;
5644  tl.push_back(model::term_description(varname, varname, true));
5645  return md.add_brick(pbr, model::varnamelist(1, varname),
5646  model::varnamelist(1, dataexpr), tl,
5647  model::mimlist(1, &mim), region);
5648  } else {
5649  std::string test_varname
5650  = "Test_" + sup_previous_and_dot_to_varname(varname);
5651  std::string expr = "(("+dataexpr+")*"+varname+")."+test_varname;
5652  size_type ib = add_linear_term(md, mim, expr, region, true, true,
5653  "Fourier-Robin", true);
5654  if (ib == size_type(-1))
5655  ib = add_nonlinear_term(md, mim, expr, region, false, false,
5656  "Fourier-Robin (nonlinear)");
5657  return ib;
5658  }
5659  }
5660 
5661  // ----------------------------------------------------------------------
5662  //
5663  // Constraint brick
5664  //
5665  // ----------------------------------------------------------------------
5666 
5667  struct have_private_data_brick : public virtual_brick {
5668 
5669  model_real_sparse_matrix rB;
5670  model_complex_sparse_matrix cB;
5671  model_real_plain_vector rL;
5672  model_complex_plain_vector cL;
5673  std::string nameL;
5674  };
5675 
5676  struct constraint_brick : public have_private_data_brick {
5677 
5678  virtual void real_pre_assembly_in_serial(const model &md, size_type,
5679  const model::varnamelist &vl,
5680  const model::varnamelist &dl,
5681  const model::mimlist &mims,
5682  model::real_matlist &matl,
5683  model::real_veclist &vecl,
5684  model::real_veclist &,
5685  size_type, build_version) const {
5686  if (MPI_IS_MASTER()) {
5687 
5688  GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5689  "Constraint brick has one and only one term");
5690  GMM_ASSERT1(mims.size() == 0,
5691  "Constraint brick need no mesh_im");
5692  GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() <= 1,
5693  "Wrong number of variables for constraint brick");
5694 
5695  bool penalized = (vl.size() == 1);
5696  const model_real_plain_vector *COEFF = 0;
5697 
5698  bool has_data = (nameL.compare("") != 0);
5699  if (has_data)
5700  GMM_ASSERT1(nameL.compare(dl.back()) == 0 &&
5701  md.variable_exists(nameL) && md.is_data(nameL),
5702  "Internal error");
5703  const model_real_plain_vector &
5704  rrL = has_data ? md.real_variable(nameL) : rL;
5705 
5706  if (penalized) {
5707  COEFF = &(md.real_variable(dl[0]));
5708  GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
5709  "Data for coefficient should be a scalar");
5710 
5711  gmm::mult(gmm::transposed(rB),
5712  gmm::scaled(rrL, gmm::abs((*COEFF)[0])), vecl[0]);
5713  gmm::mult(gmm::transposed(rB),
5714  gmm::scaled(rB, gmm::abs((*COEFF)[0])), matl[0]);
5715  } else {
5716  gmm::copy(rrL, vecl[0]);
5717  gmm::copy(rB, matl[0]);
5718  }
5719  }
5720  }
5721 
5722  virtual void complex_pre_assembly_in_serial(const model &md, size_type,
5723  const model::varnamelist &vl,
5724  const model::varnamelist &dl,
5725  const model::mimlist &mims,
5726  model::complex_matlist &matl,
5727  model::complex_veclist &vecl,
5728  model::complex_veclist &,
5729  size_type,
5730  build_version) const {
5731  if (MPI_IS_MASTER()) {
5732 
5733  GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5734  "Constraint brick has one and only one term");
5735  GMM_ASSERT1(mims.size() == 0,
5736  "Constraint brick need no mesh_im");
5737  GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() <= 1,
5738  "Wrong number of variables for constraint brick");
5739 
5740  bool penalized = (vl.size() == 1);
5741  const model_complex_plain_vector *COEFF = 0;
5742 
5743  bool has_data = (nameL.compare("") != 0);
5744  if (has_data)
5745  GMM_ASSERT1(nameL.compare(dl.back()) == 0 &&
5746  md.variable_exists(nameL) && md.is_data(nameL),
5747  "Internal error");
5748  const model_complex_plain_vector &
5749  ccL = has_data ? md.complex_variable(nameL) : cL;
5750 
5751  if (penalized) {
5752  COEFF = &(md.complex_variable(dl[0]));
5753  GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
5754  "Data for coefficient should be a scalar");
5755 
5756  gmm::mult(gmm::transposed(cB),
5757  gmm::scaled(ccL, gmm::abs((*COEFF)[0])), vecl[0]);
5758  gmm::mult(gmm::transposed(cB),
5759  gmm::scaled(cB, gmm::abs((*COEFF)[0])), matl[0]);
5760  } else {
5761  gmm::copy(ccL, vecl[0]);
5762  gmm::copy(cB, matl[0]);
5763  }
5764  }
5765  }
5766 
5767  virtual std::string declare_volume_assembly_string
5768  (const model &, size_type, const model::varnamelist &,
5769  const model::varnamelist &) const {
5770  return std::string();
5771  }
5772 
5773  constraint_brick(bool penalized) {
5774  set_flags(penalized ? "Constraint with penalization brick"
5775  : "Constraint with multipliers brick",
5776  true /* is linear*/,
5777  true /* is symmetric */, penalized /* is coercive */,
5778  true /* is real */, true /* is complex */,
5779  false /* compute each time */);
5780  }
5781 
5782  };
5783 
5784  model_real_sparse_matrix &set_private_data_brick_real_matrix
5785  (model &md, size_type indbrick) {
5786  pbrick pbr = md.brick_pointer(indbrick);
5787  md.touch_brick(indbrick);
5788  have_private_data_brick *p = dynamic_cast<have_private_data_brick *>
5789  (const_cast<virtual_brick *>(pbr.get()));
5790  GMM_ASSERT1(p, "Wrong type of brick");
5791  return p->rB;
5792  }
5793 
5794  model_real_plain_vector &set_private_data_brick_real_rhs
5795  (model &md, size_type indbrick) {
5796  pbrick pbr = md.brick_pointer(indbrick);
5797  md.touch_brick(indbrick);
5798  have_private_data_brick *p = dynamic_cast<have_private_data_brick *>
5799  (const_cast<virtual_brick *>(pbr.get()));
5800  GMM_ASSERT1(p, "Wrong type of brick");
5801  if (p->nameL.compare("") != 0) GMM_WARNING1("Rhs already set by data name");
5802  return p->rL;
5803  }
5804 
5805  model_complex_sparse_matrix &set_private_data_brick_complex_matrix
5806  (model &md, size_type indbrick) {
5807  pbrick pbr = md.brick_pointer(indbrick);
5808  md.touch_brick(indbrick);
5809  have_private_data_brick *p = dynamic_cast<have_private_data_brick *>
5810  (const_cast<virtual_brick *>(pbr.get()));
5811  GMM_ASSERT1(p, "Wrong type of brick");
5812  return p->cB;
5813  }
5814 
5815  model_complex_plain_vector &set_private_data_brick_complex_rhs
5816  (model &md, size_type indbrick) {
5817  pbrick pbr = md.brick_pointer(indbrick);
5818  md.touch_brick(indbrick);
5819  have_private_data_brick *p = dynamic_cast<have_private_data_brick *>
5820  (const_cast<virtual_brick *>(pbr.get()));
5821  GMM_ASSERT1(p, "Wrong type of brick");
5822  if (p->nameL.compare("") != 0) GMM_WARNING1("Rhs already set by data name");
5823  return p->cL;
5824  }
5825 
5826  void set_private_data_rhs
5827  (model &md, size_type indbrick, const std::string &varname) {
5828  pbrick pbr = md.brick_pointer(indbrick);
5829  md.touch_brick(indbrick);
5830  have_private_data_brick *p = dynamic_cast<have_private_data_brick *>
5831  (const_cast<virtual_brick *>(pbr.get()));
5832  GMM_ASSERT1(p, "Wrong type of brick");
5833  if (p->nameL.compare(varname) != 0) {
5834  model::varnamelist dl = md.datanamelist_of_brick(indbrick);
5835  if (p->nameL.compare("") == 0) dl.push_back(varname);
5836  else dl.back() = varname;
5837  md.change_data_of_brick(indbrick, dl);
5838  p->nameL = varname;
5839  }
5840  }
5841 
5842  size_type add_constraint_with_penalization
5843  (model &md, const std::string &varname, scalar_type penalisation_coeff) {
5844  std::string coeffname = md.new_name("penalization_on_" + varname);
5845  md.add_fixed_size_data(coeffname, 1);
5846  if (md.is_complex())
5847  md.set_complex_variable(coeffname)[0] = penalisation_coeff;
5848  else
5849  md.set_real_variable(coeffname)[0] = penalisation_coeff;
5850  pbrick pbr = std::make_shared<constraint_brick>(true);
5851  model::termlist tl;
5852  tl.push_back(model::term_description(varname, varname, true));
5853  model::varnamelist vl(1, varname);
5854  model::varnamelist dl(1, coeffname);
5855  return md.add_brick(pbr, vl, dl, tl, model::mimlist(), size_type(-1));
5856  }
5857 
5858  size_type add_constraint_with_multipliers
5859  (model &md, const std::string &varname, const std::string &multname) {
5860  pbrick pbr = std::make_shared<constraint_brick>(false);
5861  model::termlist tl;
5862  tl.push_back(model::term_description(multname, varname, true));
5863  model::varnamelist vl(1, varname);
5864  vl.push_back(multname);
5865  model::varnamelist dl;
5866  return md.add_brick(pbr, vl, dl, tl, model::mimlist(), size_type(-1));
5867  }
5868 
5869 
5870  // ----------------------------------------------------------------------
5871  //
5872  // Explicit matrix brick
5873  //
5874  // ----------------------------------------------------------------------
5875 
5876  struct explicit_matrix_brick : public have_private_data_brick {
5877 
5878  virtual void real_pre_assembly_in_serial(const model &, size_type,
5879  const model::varnamelist &vl,
5880  const model::varnamelist &dl,
5881  const model::mimlist &mims,
5882  model::real_matlist &matl,
5883  model::real_veclist &vecl,
5884  model::real_veclist &,
5885  size_type, build_version) const {
5886  GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5887  "Explicit matrix has one and only one term");
5888  GMM_ASSERT1(mims.size() == 0, "Explicit matrix need no mesh_im");
5889  GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() == 0,
5890  "Wrong number of variables for explicit matrix brick");
5891  GMM_ASSERT1(gmm::mat_ncols(rB) == gmm::mat_ncols(matl[0]) &&
5892  gmm::mat_nrows(rB) == gmm::mat_nrows(matl[0]),
5893  "Explicit matrix brick dimension mismatch ("<<
5894  gmm::mat_ncols(rB)<<"x"<<gmm::mat_nrows(rB)<<") != ("<<
5895  gmm::mat_ncols(matl[0])<<"x"<<gmm::mat_nrows(matl[0])<<")");
5896  gmm::copy(rB, matl[0]);
5897  }
5898 
5899  virtual void complex_pre_assembly_in_serial(const model &, size_type,
5900  const model::varnamelist &vl,
5901  const model::varnamelist &dl,
5902  const model::mimlist &mims,
5903  model::complex_matlist &matl,
5904  model::complex_veclist &vecl,
5905  model::complex_veclist &,
5906  size_type,
5907  build_version) const {
5908  GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5909  "Explicit matrix has one and only one term");
5910  GMM_ASSERT1(mims.size() == 0, "Explicit matrix need no mesh_im");
5911  GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() == 0,
5912  "Wrong number of variables for explicit matrix brick");
5913  gmm::copy(cB, matl[0]);
5914  }
5915 
5916  virtual std::string declare_volume_assembly_string
5917  (const model &, size_type, const model::varnamelist &,
5918  const model::varnamelist &) const {
5919  return std::string();
5920  }
5921 
5922  explicit_matrix_brick(bool symmetric_, bool coercive_) {
5923  set_flags("Explicit matrix brick",
5924  true /* is linear*/,
5925  symmetric_ /* is symmetric */, coercive_ /* is coercive */,
5926  true /* is real */, true /* is complex */,
5927  true /* is to be computed each time */);
5928  }
5929  };
5930 
5931  size_type add_explicit_matrix
5932  (model &md, const std::string &varname1, const std::string &varname2,
5933  bool issymmetric, bool iscoercive) {
5934  pbrick pbr = std::make_shared<explicit_matrix_brick>(issymmetric,
5935  iscoercive);
5936  model::termlist tl;
5937  tl.push_back(model::term_description(varname1, varname2, issymmetric));
5938  model::varnamelist vl(1, varname1);
5939  vl.push_back(varname2);
5940  model::varnamelist dl;
5941  return md.add_brick(pbr, vl, dl, tl, model::mimlist(), size_type(-1));
5942  }
5943 
5944  // ----------------------------------------------------------------------
5945  //
5946  // Explicit rhs brick
5947  //
5948  // ----------------------------------------------------------------------
5949 
5950  struct explicit_rhs_brick : public have_private_data_brick {
5951 
5952  virtual void real_pre_assembly_in_serial(const model &, size_type,
5953  const model::varnamelist &vl,
5954  const model::varnamelist &dl,
5955  const model::mimlist &mims,
5956  model::real_matlist &matl,
5957  model::real_veclist &vecl,
5958  model::real_veclist &,
5959  size_type, build_version) const {
5960  if (MPI_IS_MASTER()) {
5961  GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5962  "Explicit rhs has one and only one term");
5963  GMM_ASSERT1(mims.size() == 0, "Explicit rhs need no mesh_im");
5964  GMM_ASSERT1(vl.size() == 1 && dl.size() == 0,
5965  "Wrong number of variables for explicit rhs brick");
5966  gmm::copy(rL, vecl[0]);
5967  }
5968  }
5969 
5970  virtual void complex_pre_assembly_in_serial(const model &, size_type,
5971  const model::varnamelist &vl,
5972  const model::varnamelist &dl,
5973  const model::mimlist &mims,
5974  model::complex_matlist &matl,
5975  model::complex_veclist &vecl,
5976  model::complex_veclist &,
5977  size_type,
5978  build_version) const {
5979  if (MPI_IS_MASTER()) {
5980  GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
5981  "Explicit rhs has one and only one term");
5982  GMM_ASSERT1(mims.size() == 0, "Explicit rhs need no mesh_im");
5983  GMM_ASSERT1(vl.size() == 1 && dl.size() == 0,
5984  "Wrong number of variables for explicit rhs brick");
5985  gmm::copy(cL, vecl[0]);
5986  }
5987 
5988  }
5989 
5990  virtual std::string declare_volume_assembly_string
5991  (const model &, size_type, const model::varnamelist &,
5992  const model::varnamelist &) const {
5993  return std::string();
5994  }
5995 
5996  explicit_rhs_brick() {
5997  set_flags("Explicit rhs brick",
5998  true /* is linear*/,
5999  true /* is symmetric */, true /* is coercive */,
6000  true /* is real */, true /* is complex */,
6001  true /* is to be computed each time */);
6002  }
6003 
6004  };
6005 
6006  size_type add_explicit_rhs
6007  (model &md, const std::string &varname) {
6008  pbrick pbr = std::make_shared<explicit_rhs_brick>();
6009  model::termlist tl;
6010  tl.push_back(model::term_description(varname));
6011  model::varnamelist vl(1, varname);
6012  model::varnamelist dl;
6013  return md.add_brick(pbr, vl, dl, tl, model::mimlist(), size_type(-1));
6014  }
6015 
6016 
6017  // ----------------------------------------------------------------------
6018  //
6019  // Isotropic linearized elasticity brick
6020  //
6021  // ----------------------------------------------------------------------
6022 
6023  struct iso_lin_elasticity_new_brick : public virtual_brick {
6024 
6025  std::string expr, dataname3;
6026 
6027  void asm_real_tangent_terms(const model &md, size_type ib,
6028  const model::varnamelist &vl,
6029  const model::varnamelist &dl,
6030  const model::mimlist &mims,
6031  model::real_matlist &matl,
6032  model::real_veclist &vecl,
6033  model::real_veclist &,
6034  size_type region,
6035  build_version version) const override {
6036  GMM_ASSERT1(vl.size() == 1, "Linearized isotropic elasticity brick "
6037  "has one and only one variable");
6038  GMM_ASSERT1(matl.size() == 1, "Linearized isotropic elasticity brick "
6039  "has one and only one term");
6040  GMM_ASSERT1(mims.size() == 1, "Linearized isotropic elasticity brick "
6041  "needs one and only one mesh_im");
6042 
6043  bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0);
6044  for (size_type i = 0; i < dl.size(); ++i) {
6045  recompute_matrix = recompute_matrix ||
6046  md.is_var_newer_than_brick(dl[i], ib);
6047  }
6048 
6049  if (recompute_matrix) {
6050  // reenables disabled variables
6051  ga_workspace workspace(md, ga_workspace::inherit::ALL);
6052  workspace.add_expression(expr, *(mims[0]), region);
6053  GMM_TRACE2(name << ": generic matrix assembly");
6054  workspace.assembly(2);
6055  scalar_type alpha = scalar_type(1)
6056  / (workspace.factor_of_variable(vl[0]));
6057  const auto &R=workspace.assembled_matrix();
6058  gmm::sub_interval I = workspace.interval_of_variable(vl[0]);
6059  gmm::copy(gmm::scaled(gmm::sub_matrix(R, I, I), alpha),
6060  matl[0]);
6061  }
6062 
6063  if (dataname3.size()) { // Pre-constraints given by an "initial"
6064  // displacement u0. Means that the computed displacement will be u - u0
6065  // The displacement u0 should be discribed on the same fem as the
6066  // variable.
6067  gmm::clear(vecl[0]);
6068  gmm::mult(matl[0],
6069  gmm::scaled(md.real_variable(dataname3), scalar_type(-1)),
6070  vecl[0]);
6071  }
6072 
6073  }
6074 
6075  void real_post_assembly_in_serial(const model &md, size_type ib,
6076  const model::varnamelist &/* vl */,
6077  const model::varnamelist &/* dl */,
6078  const model::mimlist &/* mims */,
6079  model::real_matlist &/*matl*/,
6080  model::real_veclist &vecl,
6081  model::real_veclist &,
6082  size_type /*region*/,
6083  build_version) const override {
6084  md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
6085  }
6086 
6087 
6088  virtual std::string declare_volume_assembly_string
6089  (const model &, size_type, const model::varnamelist &,
6090  const model::varnamelist &) const {
6091  return expr;
6092  }
6093 
6094  iso_lin_elasticity_new_brick(const std::string &expr_,
6095  const std::string &dataname3_) {
6096  expr = expr_; dataname3 = dataname3_;
6097  set_flags("Linearized isotropic elasticity", true /* is linear*/,
6098  true /* is symmetric */, true /* is coercive */,
6099  true /* is real */, false /* is complex */);
6100  }
6101 
6102  };
6103 
6104 
6106  (model &md, const mesh_im &mim, const std::string &varname,
6107  const std::string &dataexpr1, const std::string &dataexpr2,
6108  size_type region, const std::string &dataname3) {
6109  std::string test_varname
6110  = "Test_" + sup_previous_and_dot_to_varname(varname);
6111 
6112  std::string expr1 = "((("+dataexpr1+")*(Div_"+varname+"-Div_"+dataname3
6113  +"))*Id(meshdim)+(2*("+dataexpr2+"))*(Sym(Grad_"+varname
6114  +")-Sym(Grad_"+dataname3+"))):Grad_" +test_varname;
6115  std::string expr2 = "(Div_"+varname+"*(("+dataexpr1+")*Id(meshdim))"
6116  +"+(2*("+dataexpr2+"))*Sym(Grad_"+varname+")):Grad_"+test_varname;
6117 
6118  bool is_lin;
6119  model::varnamelist vl, dl;
6120  { // reenables disabled variables
6121  ga_workspace workspace(md, ga_workspace::inherit::ALL);
6122  workspace.add_expression(expr2, mim, region);
6123  model::varnamelist vl_test1, vl_test2;
6124  is_lin = workspace.used_variables(vl, vl_test1, vl_test2, dl, 2);
6125  }
6126  if (is_lin) {
6127  pbrick pbr = std::make_shared<iso_lin_elasticity_new_brick>
6128  (expr2, dataname3);
6129  model::termlist tl;
6130  tl.push_back(model::term_description(varname,
6131  sup_previous_and_dot_to_varname(varname), true));
6132  if (dataname3.size()) dl.push_back(dataname3);
6133  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
6134  } else {
6135  return add_nonlinear_generic_assembly_brick
6136  (md, mim, dataname3.size() ? expr1 : expr2, region, false, false,
6137  "Linearized isotropic elasticity (with nonlinear dependance)");
6138  }
6139  }
6140 
6142  (model &md, const mesh_im &mim, const std::string &varname,
6143  const std::string &data_E, const std::string &data_nu,
6144  size_type region) {
6145  std::string test_varname
6146  = "Test_" + sup_previous_and_dot_to_varname(varname);
6147 
6148  std::string mu = "(("+data_E+")/(2*(1+("+data_nu+"))))";
6149  std::string lambda = "(("+data_E+")*("+data_nu+")/((1+("+data_nu+"))*(1-2*("
6150  +data_nu+"))))";
6151  std::string expr = lambda+"*Div_"+varname+"*Div_"+test_varname
6152  + "+"+mu+"*(Grad_"+varname+"+Grad_"+varname+"'):Grad_"+test_varname;
6153 
6154  bool is_lin;
6155  { // reenables disabled variables
6156  ga_workspace workspace(md, ga_workspace::inherit::ALL);
6157  workspace.add_expression(expr, mim, region);
6158  is_lin = workspace.is_linear(2);
6159  }
6160  if (is_lin) {
6161  return add_linear_term(md, mim, expr, region, false, false,
6162  "Linearized isotropic elasticity");
6163  } else {
6164  return add_nonlinear_term
6165  (md, mim, expr, region, false, false,
6166  "Linearized isotropic elasticity (with nonlinear dependance)");
6167  }
6168  }
6169 
6171  (model &md, const mesh_im &mim, const std::string &varname,
6172  const std::string &data_E, const std::string &data_nu,
6173  size_type region) {
6174  std::string test_varname
6175  = "Test_" + sup_previous_and_dot_to_varname(varname);
6176 
6177  const mesh_fem *mfu = md.pmesh_fem_of_variable(varname);
6178  GMM_ASSERT1(mfu, "The variable should be a fem variable");
6179  size_type N = mfu->linked_mesh().dim();
6180 
6181  std::string mu = "(("+data_E+")/(2*(1+("+data_nu+"))))";
6182  std::string lambda = "(("+data_E+")*("+data_nu+")/((1+("+data_nu+"))*(1-2*("
6183  +data_nu+"))))";
6184  if (N == 2)
6185  lambda = "(("+data_E+")*("+data_nu+")/((1-sqr("+data_nu+"))))";
6186  std::string expr = lambda+"*Div_"+varname+"*Div_"+test_varname
6187  + "+"+mu+"*(Grad_"+varname+"+Grad_"+varname+"'):Grad_"+test_varname;
6188 
6189  bool is_lin;
6190  { // reenables disabled variables
6191  ga_workspace workspace(md, ga_workspace::inherit::ALL);
6192  workspace.add_expression(expr, mim, region);
6193  is_lin = workspace.is_linear(2);
6194  }
6195  if (is_lin) {
6196  return add_linear_term(md, mim, expr, region, false, false,
6197  "Linearized isotropic elasticity");
6198  } else {
6199  return add_nonlinear_term
6200  (md, mim, expr, region, false, false,
6201  "Linearized isotropic elasticity (with nonlinear dependance)");
6202  }
6203  }
6204 
6205  // Tresca to be implemented with generic interpolation
6206  void compute_isotropic_linearized_Von_Mises_or_Tresca
6207  (model &md, const std::string &varname, const std::string &data_lambda,
6208  const std::string &data_mu, const mesh_fem &mf_vm,
6209  model_real_plain_vector &VM, bool tresca) {
6210 
6211  if (tresca) {
6212  const mesh_fem &mf_u = md.mesh_fem_of_variable(varname);
6213  const mesh_fem *mf_lambda = md.pmesh_fem_of_variable(data_lambda);
6214  const model_real_plain_vector *lambda=&(md.real_variable(data_lambda));
6215  const mesh_fem *mf_mu = md.pmesh_fem_of_variable(data_mu);
6216  const model_real_plain_vector *mu = &(md.real_variable(data_mu));
6217 
6218  size_type sl = gmm::vect_size(*lambda);
6219  if (mf_lambda) sl = sl * mf_lambda->get_qdim() / mf_lambda->nb_dof();
6220  size_type sm = gmm::vect_size(*mu);
6221  if (mf_mu) sm = sm * mf_mu->get_qdim() / mf_mu->nb_dof();
6222 
6223  GMM_ASSERT1(sl == 1 && sm == 1, "Bad format for Lame coefficients");
6224  GMM_ASSERT1(mf_lambda == mf_mu,
6225  "The two Lame coefficients should be described on the same "
6226  "finite element method.");
6227 
6228  if (mf_lambda) {
6230  md.real_variable(varname), VM,
6231  *mf_lambda, *lambda,
6232  *mf_lambda, *mu,
6233  tresca);
6234  } else {
6235  mf_lambda = &(classical_mesh_fem(mf_u.linked_mesh(), 0));
6236  model_real_plain_vector LAMBDA(mf_lambda->nb_dof(), (*lambda)[0]);
6237  model_real_plain_vector MU(mf_lambda->nb_dof(), (*mu)[0]);
6239  md.real_variable(varname), VM,
6240  *mf_lambda, LAMBDA,
6241  *mf_lambda, MU,
6242  tresca);
6243  }
6244  } else {
6245  // The Lambda part is not necessary for Von Mises stress ...
6246  // std::string sigma = "("+data_lambda+")*Div_"+varname+"*Id(meshdim)+("
6247  // + data_mu+")*(Grad_"+varname+"+Grad_"+varname+"')";
6248  std::string sigma_d="("+data_mu+")*(Grad_"+varname+"+Grad_"+varname+"')";
6249  std::string expr = "sqrt(3/2)*Norm(Deviator("+sigma_d+"))";
6250  ga_interpolation_Lagrange_fem(md, expr, mf_vm, VM);
6251  }
6252  }
6253 
6254 
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) {
6259  // The Lambda part is not necessary for Von Mises stress ...
6260  // std::string lambda = "(("+data_E+")*("+data_nu+")/((1+("+data_nu+"))*(1-2*("
6261  // +data_nu+"))))";
6262  std::string mu = "(("+data_E+")/(2*(1+("+data_nu+"))))";
6263  // std::string sigma = lambda+"*Div_"+varname+"*Id(meshdim)+"
6264  // + mu+"*(Grad_"+varname+"+Grad_"+varname+"')";
6265  std::string sigma_d = mu+"*(Grad_"+varname+"+Grad_"+varname+"')";
6266  std::string expr = "sqrt(3/2)*Norm(Deviator("+sigma_d+"))";
6267  ga_interpolation_Lagrange_fem(md, expr, mf_vm, VM);
6268  }
6269 
6271  (model &md, const std::string &varname, const std::string &data_E,
6272  const std::string &data_nu, const mesh_fem &mf_vm,
6273  model_real_plain_vector &VM) {
6274  // The Lambda part is not necessary for Von Mises stress ...
6275  // const mesh_fem *mfu = md.pmesh_fem_of_variable(varname);
6276  // GMM_ASSERT1(mfu, "The variable should be a fem variable");
6277  // size_type N = mfu->linked_mesh().dim();
6278  // std::string lambda = "(("+data_E+")*("+data_nu+")/((1+("+data_nu
6279  // +"))*(1-2*("+data_nu+"))))";
6280  // if (N == 2)
6281  // lambda = "(("+data_E+")*("+data_nu+")/((1-sqr("+data_nu+"))))";
6282  std::string mu = "(("+data_E+")/(2*(1+("+data_nu+"))))";
6283  // std::string sigma = lambda+"*Div_"+varname+"*Id(meshdim)+"
6284  // + mu+"*(Grad_"+varname+"+Grad_"+varname+"')";
6285  std::string sigma_d = mu+"*(Grad_"+varname+"+Grad_"+varname+"')";
6286  std::string expr = "sqrt(3/2)*Norm(Deviator("+sigma_d+"))";
6287  ga_interpolation_Lagrange_fem(md, expr, mf_vm, VM);
6288  }
6289 
6290 
6291  // --------------------------------------------------------------------
6292  //
6293  // linearized incompressibility brick (div u = 0)
6294  //
6295  // ----------------------------------------------------------------------
6296 
6297  struct linear_incompressibility_brick : public virtual_brick {
6298 
6299  virtual void asm_real_tangent_terms(const model &md, size_type /*ib*/,
6300  const model::varnamelist &vl,
6301  const model::varnamelist &dl,
6302  const model::mimlist &mims,
6303  model::real_matlist &matl,
6304  model::real_veclist &,
6305  model::real_veclist &,
6306  size_type region,
6307  build_version) const {
6308 
6309  GMM_ASSERT1((matl.size() == 1 && dl.size() == 0)
6310  || (matl.size() == 2 && dl.size() == 1),
6311  "Wrong term and/or data number for Linear incompressibility "
6312  "brick.");
6313  GMM_ASSERT1(mims.size() == 1, "Linear incompressibility brick need one "
6314  "and only one mesh_im");
6315  GMM_ASSERT1(vl.size() == 2, "Wrong number of variables for linear "
6316  "incompressibility brick");
6317 
6318  bool penalized = (dl.size() == 1);
6319  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
6320  const mesh_fem &mf_p = md.mesh_fem_of_variable(vl[1]);
6321  const mesh_im &mim = *mims[0];
6322  const model_real_plain_vector *COEFF = 0;
6323  const mesh_fem *mf_data = 0;
6324 
6325  if (penalized) {
6326  COEFF = &(md.real_variable(dl[0]));
6327  mf_data = md.pmesh_fem_of_variable(dl[0]);
6328  size_type s = gmm::vect_size(*COEFF);
6329  if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
6330  GMM_ASSERT1(s == 1, "Bad format for the penalization parameter");
6331  }
6332 
6333  mesh_region rg(region);
6334  mim.linked_mesh().intersect_with_mpi_region(rg);
6335 
6336  GMM_TRACE2("Stokes term assembly");
6337  gmm::clear(matl[0]);
6338  asm_stokes_B(matl[0], mim, mf_u, mf_p, rg);
6339 
6340  if (penalized) {
6341  gmm::clear(matl[1]);
6342  if (mf_data) {
6343  asm_mass_matrix_param(matl[1], mim, mf_p, *mf_data, *COEFF, rg);
6344  gmm::scale(matl[1], scalar_type(-1));
6345  }
6346  else {
6347  asm_mass_matrix(matl[1], mim, mf_p, rg);
6348  gmm::scale(matl[1], -(*COEFF)[0]);
6349  }
6350  }
6351 
6352  }
6353 
6354 
6355  virtual void real_post_assembly_in_serial(const model &, size_type,
6356  const model::varnamelist &,
6357  const model::varnamelist &/*dl*/,
6358  const model::mimlist &/*mims*/,
6359  model::real_matlist &/*matl*/,
6360  model::real_veclist &,
6361  model::real_veclist &,
6362  size_type /*region*/,
6363  build_version) const
6364  { }
6365 
6366 
6367  linear_incompressibility_brick() {
6368  set_flags("Linear incompressibility brick",
6369  true /* is linear*/,
6370  true /* is symmetric */, false /* is coercive */,
6371  true /* is real */, false /* is complex */);
6372  }
6373 
6374  };
6375 
6377  (model &md, const mesh_im &mim, const std::string &varname,
6378  const std::string &multname, size_type region,
6379  const std::string &dataexpr) {
6380 #if 0
6381  pbrick pbr = std::make_shared<linear_incompressibility_brick>();
6382  model::termlist tl;
6383  tl.push_back(model::term_description(multname, varname, true));
6384  model::varnamelist vl(1, varname);
6385  vl.push_back(multname);
6386  model::varnamelist dl;
6387  if (dataname.size()) {
6388  dl.push_back(dataexpr);
6389  tl.push_back(model::term_description(multname, multname, true));
6390  }
6391  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
6392 #else
6393  std::string test_varname
6394  = "Test_" + sup_previous_and_dot_to_varname(varname);
6395  std::string test_multname
6396  = "Test_" + sup_previous_and_dot_to_varname(multname);
6397  std::string expr;
6398  if (dataexpr.size())
6399  expr = "-"+multname+"*Div_"+test_varname + "-"+test_multname
6400  +"*Div_"+varname+"+(("+dataexpr+")*"+multname+")*"+test_multname;
6401  else
6402  expr = "-"+multname+"*Div_"+test_varname + "-"+test_multname
6403  +"*Div_"+varname;
6404  size_type ib = add_linear_term(md, mim, expr, region, true, true,
6405  "Linear incompressibility", true);
6406  if (ib == size_type(-1))
6407  ib = add_nonlinear_term
6408  (md, mim, expr, region, false, false,
6409  "Linear incompressibility (with nonlinear dependance)");
6410  return ib;
6411 #endif
6412  }
6413 
6414 
6415 
6416  // ----------------------------------------------------------------------
6417  //
6418  // Mass brick
6419  //
6420  // ----------------------------------------------------------------------
6421 
6422  struct mass_brick : public virtual_brick {
6423 
6424  virtual void asm_real_tangent_terms(const model &md, size_type,
6425  const model::varnamelist &vl,
6426  const model::varnamelist &dl,
6427  const model::mimlist &mims,
6428  model::real_matlist &matl,
6429  model::real_veclist &,
6430  model::real_veclist &,
6431  size_type region,