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