GetFEM  5.5
getfem_generic_assembly.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2013-2026 Yves Renard
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program. If not, see https://www.gnu.org/licenses/.
19 
20  As a special exception, you may use this file as it is a part of a free
21  software library without restriction. Specifically, if other files
22  instantiate templates or use macros or inline functions from this file,
23  or you compile this file and link it with other files to produce an
24  executable, this file does not by itself cause the resulting executable
25  to be covered by the GNU Lesser General Public License. This exception
26  does not however invalidate any other reasons why the executable file
27  might be covered by the GNU Lesser General Public License.
28 
29 ===========================================================================*/
30 
31 /** @file getfem_generic_assembly.h
32  @author Yves Renard <Yves.Renard@insa-lyon.fr>
33  @date November 18, 2013.
34  @brief A language for generic assembly of pde boundary value problems.
35  */
36 
37 
38 #ifndef GETFEM_GENERIC_ASSEMBLY_H__
39 #define GETFEM_GENERIC_ASSEMBLY_H__
40 
41 #include <map>
42 #include "getfem/getfem_interpolation.h"
43 #include "getfem/getfem_mesh_slice.h"
44 
45 
46 #ifdef _WIN32
47 #include <limits>
48 #if defined(INFINITY)
49 #undef INFINITY
50 #endif
51 #define INFINITY std::numeric_limits<scalar_type>::infinity()
52 #endif
53 
54 namespace getfem {
55 
56  struct ga_tree;
57  class model;
58  class ga_workspace;
59 
60  typedef gmm::rsvector<scalar_type> model_real_sparse_vector;
61  typedef gmm::rsvector<complex_type> model_complex_sparse_vector;
62  typedef std::vector<scalar_type> model_real_plain_vector;
63  typedef std::vector<complex_type> model_complex_plain_vector;
64 
65  typedef gmm::col_matrix<model_real_sparse_vector> model_real_sparse_matrix;
66  typedef gmm::col_matrix<model_complex_sparse_vector>
67  model_complex_sparse_matrix;
68 
69  typedef gmm::row_matrix<model_real_sparse_vector>
70  model_real_row_sparse_matrix;
71  typedef gmm::row_matrix<model_complex_sparse_vector>
72  model_complex_row_sparse_matrix;
73 
74  // 0 : ok
75  // 1 : function or operator name or "X"
76  // 2 : reserved prefix Grad, Hess, Div, Test and Test2
77  int ga_check_name_validity(const std::string &name);
78 
79  //=========================================================================
80  // Virtual interpolate_transformation object.
81  //=========================================================================
82 
83  struct var_trans_pair {
84  std::string varname, transname;
85  bool operator <(const var_trans_pair &vt) const {
86  return (varname < vt.varname) ||
87  (!(varname > vt.varname) && transname < vt.transname);
88  }
89  var_trans_pair() : varname(), transname() {}
90  var_trans_pair(const std::string &v, const std::string &t)
91  : varname(v), transname(t) {}
92  };
93 
94  class APIDECL virtual_interpolate_transformation {
95 
96  public:
97  virtual void extract_variables
98  (const ga_workspace &workspace, std::set<var_trans_pair> &vars,
99  bool ignore_data, const mesh &m,
100  const std::string &interpolate_name) const = 0;
101  virtual void init(const ga_workspace &workspace) const = 0;
102  virtual int transform
103  (const ga_workspace &workspace, const mesh &m,
104  fem_interpolation_context &ctx_x, const base_small_vector &Normal,
105  const mesh **m_t, size_type &cv, short_type &face_num,
106  base_node &P_ref, base_small_vector &N_y,
107  std::map<var_trans_pair, base_tensor> &derivatives,
108  bool compute_derivatives) const = 0;
109  virtual void finalize() const = 0;
110  virtual std::string expression() const { return std::string(); }
111 
112  virtual ~virtual_interpolate_transformation() {}
113  };
114 
115  typedef std::shared_ptr<const virtual_interpolate_transformation>
116  pinterpolate_transformation;
117 
118  //=========================================================================
119  // Virtual elementary_transformation object.
120  //=========================================================================
121 
122  class APIDECL virtual_elementary_transformation {
123 
124  public:
125 
126  virtual void give_transformation(const mesh_fem &mf1, const mesh_fem &mf2,
127  size_type cv, base_matrix &M) const = 0;
128  virtual ~virtual_elementary_transformation() {}
129  };
130 
131  typedef std::shared_ptr<const virtual_elementary_transformation>
132  pelementary_transformation;
133 
134  //=========================================================================
135  // Virtual secondary_domain object.
136  //=========================================================================
137 
138  class APIDECL virtual_secondary_domain {
139  protected:
140  const mesh_im &mim_;
141  const mesh_region region;
142 
143  public:
144 
145  const mesh_im &mim() const { return mim_; }
146  virtual const mesh_region &give_region(const mesh &m,
147  size_type cv, short_type f) const = 0;
148  // virtual void init(const ga_workspace &workspace) const = 0;
149  // virtual void finalize() const = 0;
150 
151  virtual_secondary_domain(const mesh_im &mim__, const mesh_region &region_)
152  : mim_(mim__), region(region_) {}
153  virtual ~virtual_secondary_domain() {}
154  };
155 
156  typedef std::shared_ptr<const virtual_secondary_domain> psecondary_domain;
157 
158  //=========================================================================
159  // Structure dealing with macros.
160  //=========================================================================
161 
162  class ga_macro {
163 
164  protected:
165  ga_tree *ptree;
166  std::string macro_name_;
167  size_type nbp;
168 
169  public:
170  ga_macro();
171  ga_macro(const std::string &name, const ga_tree &t, size_type nbp_);
172  ga_macro(const ga_macro &);
173  ~ga_macro();
174  ga_macro &operator =(const ga_macro &);
175 
176  const std::string &name() const { return macro_name_; }
177  std::string &name() { return macro_name_; }
178  size_type nb_params() const { return nbp; }
179  size_type &nb_params() { return nbp; }
180  const ga_tree& tree() const { return *ptree; }
181  ga_tree& tree() { return *ptree; }
182  };
183 
184 
185  class ga_macro_dictionary {
186 
187  protected:
188  const ga_macro_dictionary *parent;
189  std::map<std::string, ga_macro> macros;
190 
191  public:
192  bool macro_exists(const std::string &name) const;
193  const ga_macro &get_macro(const std::string &name) const;
194 
195  void add_macro(const ga_macro &gam);
196  void add_macro(const std::string &name, const std::string &expr);
197  void del_macro(const std::string &name);
198 
199  ga_macro_dictionary() : parent(0) {}
200  ga_macro_dictionary(bool, const ga_macro_dictionary& gamd)
201  : parent(&gamd) {}
202 
203  };
204 
205  //=========================================================================
206  // Structure dealing with predefined operators.
207  //=========================================================================
208 
209 
210  struct ga_nonlinear_operator {
211 
212  typedef std::vector<const base_tensor *> arg_list;
213 
214  virtual bool result_size(const arg_list &args,
215  bgeot::multi_index &sizes) const = 0;
216 
217  virtual void value(const arg_list &args, base_tensor &result) const = 0;
218 
219  virtual void derivative(const arg_list &args, size_type i,
220  base_tensor &result) const = 0;
221 
222  virtual void second_derivative(const arg_list &args, size_type i,
223  size_type j, base_tensor &result) const = 0;
224 
225  virtual ~ga_nonlinear_operator() {}
226  };
227 
228  struct ga_predef_operator_tab {
229  typedef std::map<std::string, std::shared_ptr<ga_nonlinear_operator>> T;
230  T tab;
231 
232  void add_method(const std::string &name,
233  const std::shared_ptr<ga_nonlinear_operator> &pt)
234  { tab[name] = pt; }
235  ga_predef_operator_tab();
236  };
237 
238  //=========================================================================
239  // For user predefined scalar functions.
240  //=========================================================================
241 
242  typedef scalar_type (*pscalar_func_onearg)(scalar_type);
243  typedef scalar_type (*pscalar_func_twoargs)(scalar_type, scalar_type);
244 
245  void ga_define_function(const std::string &name, size_type nb_args,
246  const std::string &expr, const std::string &der1="",
247  const std::string &der2="");
248  void ga_define_function(const std::string &name, pscalar_func_onearg f,
249  const std::string &der1="");
250  void ga_define_function(const std::string &name, pscalar_func_twoargs f2,
251  const std::string &der1="",
252  const std::string &der2="");
253 
254  void ga_undefine_function(const std::string &name);
255  bool ga_function_exists(const std::string &name);
256 
257  //=========================================================================
258  // Structure dealing with user defined environment : constant, variables,
259  // functions, operators.
260  //=========================================================================
261 
262  class ga_workspace {
263 
264  const model *md;
265  const ga_workspace *parent_workspace;
266  bool with_parent_variables;
267  size_type nb_prim_dof, nb_intern_dof, first_intern_dof, nb_tmp_dof;
268 
269  void init();
270 
271  struct var_description {
272 
273  bool is_variable;
274  const mesh_fem *mf;
275  const im_data *imd;
276  gmm::sub_interval I;
277  const model_real_plain_vector *V;
278  bgeot::multi_index qdims; // For data having a qdim different than the
279  // qdim of the fem or im_data (dim per dof for
280  // dof data) and for constant variables.
281  bool is_internal;
282 
283  size_type qdim() const {
284  size_type q = 1;
285  for (size_type i = 0; i < qdims.size(); ++i) q *= qdims[i];
286  return q;
287  }
288 
289  var_description(bool is_var, const mesh_fem *mf_, const im_data *imd_,
290  gmm::sub_interval I_, const model_real_plain_vector *V_,
291  size_type Q, bool is_intern_=false)
292  : is_variable(is_var), mf(mf_), imd(imd_),
293  I(I_), V(V_), qdims(1), is_internal(is_intern_)
294  {
295  GMM_ASSERT1(Q > 0, "Bad dimension");
296  qdims[0] = Q;
297  }
298  };
299 
300  public:
301 
302  enum operation_type {ASSEMBLY,
303  PRE_ASSIGNMENT,
304  POST_ASSIGNMENT};
305 
306  struct tree_description { // CAUTION: Specific copy constructor
307  size_type order; // 0 : potential
308  // 1 : residual
309  // 2 : tangent operator
310  // -1 : any
311  operation_type operation;
312  std::string varname_interpolation; // Where to interpolate
313  std::string name_test1, name_test2;
314  std::string interpolate_name_test1, interpolate_name_test2;
315  std::string secondary_domain;
316  const mesh_im *mim;
317  const mesh *m;
318  const mesh_region *rg;
319  ga_tree *ptree;
320  tree_description()
321  : operation(ASSEMBLY), varname_interpolation(""),
322  name_test1(""), name_test2(""),
323  interpolate_name_test1(""), interpolate_name_test2(""),
324  mim(0), m(0), rg(0), ptree(0) {}
325  void copy(const tree_description& td);
326  tree_description(const tree_description& td) { copy(td); }
327  tree_description &operator =(const tree_description& td);
328  ~tree_description();
329  };
330 
331  mutable std::set<var_trans_pair> test1, test2;
332  var_trans_pair selected_test1, selected_test2;
333 
334  private:
335 
336  // mesh regions
337  std::map<const mesh *, std::list<mesh_region> > registered_mesh_regions;
338 
339  const mesh_region &register_region(const mesh &m, const mesh_region &rg);
340 
341  // variables and variable groups
342  typedef std::map<std::string, var_description> VAR_SET;
343  VAR_SET variables;
344 
345  std::map<std::string, gmm::sub_interval> reenabled_var_intervals,
346  tmp_var_intervals;
347 
348  std::map<std::string, pinterpolate_transformation> transformations;
349  std::map<std::string, pelementary_transformation> elem_transformations;
350  std::map<std::string, psecondary_domain> secondary_domains;
351 
352  std::vector<tree_description> trees;
353 
354  std::map<std::string, std::vector<std::string> > variable_groups;
355 
356  ga_macro_dictionary macro_dict;
357 
358  // Adds a tree to the workspace. The argument tree is consumed by the
359  // function and cannot be reused afterwards.
360  void add_tree(ga_tree &tree, const mesh &m, const mesh_im &mim,
361  const mesh_region &rg, size_type add_derivative_order,
362  bool scalar_expr, operation_type op_type=ASSEMBLY,
363  const std::string varname_interpolation="");
364 
365  std::shared_ptr<model_real_sparse_matrix> K, KQJpr;
366  std::shared_ptr<base_vector> V; // reduced residual vector (primary vars + internal vars)
367  // after condensation it partially holds the condensed residual
368  // and the internal solution
369  model_real_sparse_matrix col_unreduced_K,
370  row_unreduced_K,
371  row_col_unreduced_K;
372  base_vector unreduced_V, cached_V;
373  base_tensor assemb_t;
374  bool include_empty_int_pts = false;
375 
376  public:
377  // setter functions
378  void set_assembled_matrix(model_real_sparse_matrix &K_) {
379  K = std::shared_ptr<model_real_sparse_matrix>
380  (std::shared_ptr<model_real_sparse_matrix>(), &K_); // alias
381  }
382  void set_assembled_vector(base_vector &V_) {
383  V = std::shared_ptr<base_vector>
384  (std::shared_ptr<base_vector>(), &V_); // alias
385  }
386  // getter functions
387  const model_real_sparse_matrix &assembled_matrix() const { return *K; }
388  model_real_sparse_matrix &assembled_matrix() { return *K; }
389  const base_vector &assembled_vector() const { return *V; }
390  base_vector &assembled_vector() { return *V; }
391  const base_vector &cached_vector() const { return cached_V; }
392  const base_tensor &assembled_tensor() const { return assemb_t; }
393  base_tensor &assembled_tensor() { return assemb_t; }
394  const scalar_type &assembled_potential() const {
395  GMM_ASSERT1(assemb_t.size() == 1, "Bad result size");
396  return assemb_t[0];
397  }
398  scalar_type &assembled_potential() {
399  GMM_ASSERT1(assemb_t.size() == 1, "Bad result size");
400  return assemb_t[0];
401  }
402  model_real_sparse_matrix &row_unreduced_matrix()
403  { return row_unreduced_K; }
404  model_real_sparse_matrix &col_unreduced_matrix()
405  { return col_unreduced_K; }
406  model_real_sparse_matrix &row_col_unreduced_matrix()
407  { return row_col_unreduced_K; }
408  base_vector &unreduced_vector() { return unreduced_V; }
409  // setter function for condensation matrix
410  void set_internal_coupling_matrix(model_real_sparse_matrix &KQJpr_) {
411  KQJpr = std::shared_ptr<model_real_sparse_matrix>
412  (std::shared_ptr<model_real_sparse_matrix>(), &KQJpr_); // alias
413  }
414  /** getter function for condensation matrix
415  Its size is (nb_primary_dof()+nb_internal_dof()) x nb_primary_dof()
416  but its first nb_primary_dof() rows are empty */
417  const model_real_sparse_matrix &internal_coupling_matrix() const
418  { return *KQJpr; }
419  model_real_sparse_matrix &internal_coupling_matrix() { return *KQJpr; }
420 
421  /** Add an expression, perform the semantic analysis, split into
422  * terms in separated test functions, derive if necessary to obtain
423  * the tangent terms. Return the maximal order found in the expression.
424  */
425  size_type add_expression(const std::string &expr, const mesh_im &mim,
426  const mesh_region &rg=mesh_region::all_convexes(),
427  size_type add_derivative_order = 2,
428  const std::string &secondary_dom = "");
429  /* Internal use */
430  void add_function_expression(const std::string &expr);
431  /* Internal use */
432  void add_interpolation_expression
433  (const std::string &expr, const mesh &m,
434  const mesh_region &rg = mesh_region::all_convexes());
435  void add_interpolation_expression
436  (const std::string &expr, const mesh_im &mim,
437  const mesh_region &rg = mesh_region::all_convexes());
438  void add_assignment_expression
439  (const std::string &dataname, const std::string &expr,
440  const mesh_region &rg_ = mesh_region::all_convexes(),
441  size_type order = 1, bool before = false);
442 
443  /** Delete all previously added expressions. */
444  void clear_expressions();
445 
446  /** Print some information about all previously added expressions. */
447  void print(std::ostream &str);
448 
449  size_type nb_trees() const { return trees.size(); }
450  const tree_description &tree_info(size_type i) const { return trees[i]; }
451 
452  // variables and variable groups
453  void add_fem_variable(const std::string &name, const mesh_fem &mf,
454  const gmm::sub_interval &I,
455  const model_real_plain_vector &VV);
456  void add_im_variable(const std::string &name, const im_data &imd,
457  const gmm::sub_interval &I,
458  const model_real_plain_vector &VV);
459  void add_internal_im_variable(const std::string &name, const im_data &imd,
460  const gmm::sub_interval &I,
461  const model_real_plain_vector &VV);
462  void add_fixed_size_variable(const std::string &name,
463  const gmm::sub_interval &I,
464  const model_real_plain_vector &VV);
465  void add_fem_constant(const std::string &name, const mesh_fem &mf,
466  const model_real_plain_vector &VV);
467  void add_fixed_size_constant(const std::string &name,
468  const model_real_plain_vector &VV);
469  void add_im_data(const std::string &name, const im_data &imd,
470  const model_real_plain_vector &VV);
471 
472  bool used_variables(std::vector<std::string> &vl,
473  std::vector<std::string> &vl_test1,
474  std::vector<std::string> &vl_test2,
475  std::vector<std::string> &dl,
476  size_type order);
477  bool is_linear(size_type order);
478 
479  bool variable_exists(const std::string &name) const;
480 
481  bool is_internal_variable(const std::string &name) const;
482 
483  const std::string &variable_in_group(const std::string &group_name,
484  const mesh &m) const;
485 
486  void define_variable_group(const std::string &group_name,
487  const std::vector<std::string> &nl);
488 
489  bool variable_group_exists(const std::string &name) const;
490 
491  bool variable_or_group_exists(const std::string &name) const
492  { return variable_exists(name) || variable_group_exists(name); }
493 
494  const std::vector<std::string> &
495  variable_group(const std::string &group_name) const;
496 
497  const std::string& first_variable_of_group(const std::string &name) const;
498 
499  bool is_constant(const std::string &name) const;
500 
501  bool is_disabled_variable(const std::string &name) const;
502 
503  const scalar_type &factor_of_variable(const std::string &name) const;
504 
505  const gmm::sub_interval &
506  interval_of_variable(const std::string &name) const;
507 
508  const mesh_fem *associated_mf(const std::string &name) const;
509 
510  const im_data *associated_im_data(const std::string &name) const;
511 
512  size_type qdim(const std::string &name) const;
513 
514  bgeot::multi_index qdims(const std::string &name) const;
515 
516  const model_real_plain_vector &value(const std::string &name) const;
517  scalar_type get_time_step() const;
518 
519  // macros
520  bool macro_exists(const std::string &name) const
521  { return macro_dict.macro_exists(name); }
522 
523  void add_macro(const std::string &name, const std::string &expr)
524  { macro_dict.add_macro(name, expr); }
525 
526  void del_macro(const std::string &name) { macro_dict.del_macro(name); }
527 
528  const std::string& get_macro(const std::string &name) const;
529 
530  const ga_macro_dictionary &macro_dictionary() const { return macro_dict; }
531 
532 
533  // interpolate and elementary transformations
534  void add_interpolate_transformation(const std::string &name,
535  pinterpolate_transformation ptrans);
536 
537  bool interpolate_transformation_exists(const std::string &name) const;
538 
539  pinterpolate_transformation
540  interpolate_transformation(const std::string &name) const;
541 
542  void add_elementary_transformation(const std::string &name,
543  pelementary_transformation ptrans)
544  { elem_transformations[name] = ptrans; }
545 
546  bool elementary_transformation_exists(const std::string &name) const;
547 
548  pelementary_transformation
549  elementary_transformation(const std::string &name) const;
550 
551  void add_secondary_domain(const std::string &name,
552  psecondary_domain psecdom);
553 
554  bool secondary_domain_exists(const std::string &name) const;
555 
556  psecondary_domain secondary_domain(const std::string &name) const;
557 
558  // extract terms
559  std::string extract_constant_term(const mesh &m);
560  std::string extract_order1_term(const std::string &varname);
561  std::string extract_order0_term();
562  std::string extract_Neumann_term(const std::string &varname);
563 
564  void assembly(size_type order, bool condensation=false);
565 
566  void set_include_empty_int_points(bool include);
567  bool include_empty_int_points() const;
568 
569  size_type nb_primary_dof() const { return nb_prim_dof; }
570  size_type nb_internal_dof() const { return nb_intern_dof; }
571  size_type first_internal_dof() const { return first_intern_dof; }
572  size_type nb_temporary_dof() const { return nb_tmp_dof; }
573 
574  void add_temporary_interval_for_unreduced_variable(const std::string &name);
575 
576  void clear_temporary_variable_intervals() {
577  tmp_var_intervals.clear();
578  nb_tmp_dof = 0;
579  }
580 
581  const gmm::sub_interval &
582  temporary_interval_of_variable(const std::string &name) const {
583  static const gmm::sub_interval empty_interval;
584  const auto it = tmp_var_intervals.find(name);
585  return (it != tmp_var_intervals.end()) ? it->second : empty_interval;
586  }
587 
588  enum class inherit { NONE, ENABLED, ALL };
589 
590  explicit ga_workspace(const getfem::model &md_,
591  const inherit var_inherit=inherit::ENABLED);
592  explicit ga_workspace(const ga_workspace &gaw, // compulsory 2nd arg to avoid
593  const inherit var_inherit); // conflict with copy constructor
594  ga_workspace();
595  ~ga_workspace();
596 
597  };
598 
599  // Small tool to make basic substitutions into an assembly string
600  std::string ga_substitute(const std::string &expr,
601  const std::map<std::string, std::string> &dict);
602 
603  inline std::string ga_substitute(const std::string &expr,
604  const std::string &o1,const std::string &s1) {
605  std::map<std::string, std::string> dict;
606  dict[o1] = s1;
607  return ga_substitute(expr, dict);
608  }
609 
610  inline std::string ga_substitute(const std::string &expr,
611  const std::string &o1,const std::string &s1,
612  const std::string &o2,const std::string &s2) {
613  std::map<std::string, std::string> dict;
614  dict[o1] = s1; dict[o2] = s2;
615  return ga_substitute(expr, dict);
616  }
617 
618  inline std::string ga_substitute(const std::string &expr,
619  const std::string &o1,const std::string &s1,
620  const std::string &o2,const std::string &s2,
621  const std::string &o3,const std::string &s3) {
622  std::map<std::string, std::string> dict;
623  dict[o1] = s1; dict[o2] = s2; dict[o3] = s3;
624  return ga_substitute(expr, dict);
625  }
626 
627  inline std::string ga_substitute(const std::string &expr,
628  const std::string &o1,const std::string &s1,
629  const std::string &o2,const std::string &s2,
630  const std::string &o3,const std::string &s3,
631  const std::string &o4,const std::string &s4) {
632  std::map<std::string, std::string> dict;
633  dict[o1] = s1; dict[o2] = s2; dict[o3] = s3; dict[o4] = s4;
634  return ga_substitute(expr, dict);
635  }
636 
637 
638  //=========================================================================
639  // Intermediate structure for user function manipulation
640  //=========================================================================
641 
642  struct ga_instruction_set;
643 
644  class ga_function {
645  mutable ga_workspace local_workspace;
646  std::string expr;
647  mutable ga_instruction_set *gis;
648 
649  public:
650  ga_function() : local_workspace(), expr(""), gis(0) {}
651  ga_function(const model &md, const std::string &e);
652  ga_function(const ga_workspace &workspace_, const std::string &e);
653  ga_function(const std::string &e);
654  ga_function(const ga_function &gaf);
655  ga_function &operator =(const ga_function &gaf);
656  ~ga_function();
657  const std::string &expression() const { return expr; }
658  const base_tensor &eval() const;
659  void derivative(const std::string &variable);
660  void compile() const;
661  ga_workspace &workspace() const { return local_workspace; }
662 
663  };
664 
665  //=========================================================================
666  // Intermediate structure for interpolation functions
667  //=========================================================================
668 
669  struct ga_interpolation_context {
670 
671  virtual bgeot::pstored_point_tab
672  ppoints_for_element(size_type cv, short_type f,
673  std::vector<size_type> &ind) const = 0;
674  inline const bgeot::stored_point_tab &points_for_element
675  (size_type cv, short_type f, std::vector<size_type> &ind) const
676  { return *ppoints_for_element(cv, f, ind); }
677  virtual bool use_pgp(size_type cv) const = 0;
678  virtual bool use_mim() const = 0;
679  virtual void store_result(size_type cv, size_type i, base_tensor &t) = 0;
680  virtual void finalize() = 0;
681  virtual const mesh &linked_mesh() = 0;
682  virtual ~ga_interpolation_context() {}
683  };
684 
685  //=========================================================================
686  // Interpolation functions
687  //=========================================================================
688 
689  void ga_interpolation(ga_workspace &workspace,
690  ga_interpolation_context &gic);
691 
692  void ga_interpolation_Lagrange_fem
693  (ga_workspace &workspace, const mesh_fem &mf, base_vector &result);
694 
695  void ga_interpolation_Lagrange_fem
696  (const getfem::model &md, const std::string &expr, const mesh_fem &mf,
697  base_vector &result, const mesh_region &rg=mesh_region::all_convexes());
698 
699  void ga_interpolation_mti
700  (const getfem::model &md, const std::string &expr, mesh_trans_inv &mti,
701  base_vector &result, const mesh_region &rg=mesh_region::all_convexes(),
702  int extrapolation = 0,
703  const mesh_region &rg_source=mesh_region::all_convexes(),
704  size_type nbdof_ = size_type(-1));
705 
706  void ga_interpolation_im_data
707  (ga_workspace &workspace, const im_data &imd, base_vector &result);
708 
709  void ga_interpolation_im_data
710  (const getfem::model &md, const std::string &expr, const im_data &imd,
711  base_vector &result, const mesh_region &rg=mesh_region::all_convexes());
712 
713  void ga_interpolation_mesh_slice
714  (ga_workspace &workspace, const stored_mesh_slice &sl, base_vector &result);
715 
716  void ga_interpolation_mesh_slice
717  (const getfem::model &md, const std::string &expr, const stored_mesh_slice &sl,
718  base_vector &result, const mesh_region &rg=mesh_region::all_convexes());
719 
720 
721  //=========================================================================
722  // Local projection functions
723  //=========================================================================
724 
725  /** Make an elementwise L2 projection of an expression with respect
726  to the mesh_fem `mf`. This mesh_fem has to be a discontinuous one.
727  The expression has to be valid according to the high-level generic
728  assembly language possibly including references to the variables
729  and data of the model.
730  */
731  void ga_local_projection(const getfem::model &md, const mesh_im &mim,
732  const std::string &expr, const mesh_fem &mf,
733  base_vector &result,
734  const mesh_region &rg=mesh_region::all_convexes());
735 
736  //=========================================================================
737  // Interpolate transformations
738  //=========================================================================
739 
740  /** Add a transformation to a workspace `workspace` or a model `md` mapping
741  point in mesh `source_mesh` to mesh `target_mesh`, optionally restricted
742  to the region `target_region`. The transformation is defined by the
743  expression `expr`, which has to be in the high-level generic assembly
744  syntax and may contain some variables of the workspace/model.
745  CAUTION: For the moment, the derivative of the transformation with
746  respect to any of these variables is not taken into account in the model
747  solve.
748  */
749  void add_interpolate_transformation_from_expression
750  (ga_workspace &workspace, const std::string &transname,
751  const mesh &source_mesh, const mesh &target_mesh, const std::string &expr);
752  void add_interpolate_transformation_from_expression
753  (ga_workspace &workspace, const std::string &transname,
754  const mesh &source_mesh, const mesh &target_mesh,
755  size_type target_region, const std::string &expr);
756  void add_interpolate_transformation_from_expression
757  (model &md, const std::string &transname,
758  const mesh &source_mesh, const mesh &target_mesh, const std::string &expr);
759  void add_interpolate_transformation_from_expression
760  (model &md, const std::string &transname,
761  const mesh &source_mesh, const mesh &target_mesh,
762  size_type target_region, const std::string &expr);
763 
764  /** Add a transformation to the workspace that creates an identity mapping
765  between two meshes in deformed state. Conceptually, it can be viewed
766  as a transformation from expression Xsource + Usource - Utarget,
767  except such an expression cannot be used directly in the transformation
768  from expression (function above), as Utarget needs to be interpolated
769  though an inversion of the transformation of the target domain.
770  Thread safe if added to thread local workspace.
771  */
772  void add_interpolate_transformation_on_deformed_domains
773  (ga_workspace &workspace, const std::string &transname,
774  const mesh &source_mesh, const std::string &source_displacements,
775  const mesh_region &source_region, const mesh &target_mesh,
776  const std::string &target_displacements, const mesh_region &target_region);
777 
778  /** The same as above, but adding transformation to the model.
779  Note, this version is not thread safe.*/
780  void add_interpolate_transformation_on_deformed_domains
781  (model &md, const std::string &transname,
782  const mesh &source_mesh, const std::string &source_displacements,
783  const mesh_region &source_region, const mesh &target_mesh,
784  const std::string &target_displacements, const mesh_region &target_region);
785 
786  /** Create a new instance of a transformation corresponding to the
787  interpolation on the neighbor element. Can only be applied to the
788  computation on some internal faces of a mesh.
789  (mainly for internal use in the constructor of getfem::model)
790  */
791  pinterpolate_transformation interpolate_transformation_neighbor_instance();
792 
793  /* Add a special interpolation transformation which represents the identity
794  transformation but allows to evaluate the expression on another element
795  than the current element by polynomial extrapolation. It is used for
796  stabilization term in fictitious domain applications. the map elt_cor
797  list the element concerned by the transformation and associate them
798  to the element on which the extrapolation has to be made. If an element
799  is not listed in elt_cor the evaluation is just made on the current
800  element.
801  */
802  void add_element_extrapolation_transformation
803  (model &md, const std::string &name, const mesh &sm,
804  std::map<size_type, size_type> &elt_corr);
805 
806  void add_element_extrapolation_transformation
807  (ga_workspace &workspace, const std::string &name, const mesh &sm,
808  std::map<size_type, size_type> &elt_corr);
809 
810  /* Change the correspondence map of an element extrapolation interpolate
811  transformation.
812  */
813  void set_element_extrapolation_correspondence
814  (model &md, const std::string &name,
815  std::map<size_type, size_type> &elt_corr);
816 
817  void set_element_extrapolation_correspondence
818  (ga_workspace &workspace, const std::string &name,
819  std::map<size_type, size_type> &elt_corr);
820 
821  //=========================================================================
822  // Secondary domains
823  //=========================================================================
824 
825  void add_standard_secondary_domain
826  (model &md, const std::string &name, const mesh_im &mim,
827  const mesh_region &rg=mesh_region::all_convexes());
828 
829  void add_standard_secondary_domain
830  (ga_workspace &workspace, const std::string &name, const mesh_im &mim,
831  const mesh_region &rg=mesh_region::all_convexes());
832 
833 
834 } /* end of namespace getfem. */
835 
836 
837 #endif /* GETFEM_GENERIC_ASSEMBLY_H__ */
getfem::mesh_region::all_convexes
static mesh_region all_convexes()
provide a default value for the mesh_region parameters of assembly procedures etc.
Definition: getfem_mesh_region.h:141
getfem::model
`‘Model’' variables store the variables, the data and the description of a model.
Definition: getfem_models.h:114
gmm::rsvector
sparse vector built upon std::vector.
Definition: gmm_vector.h:995
getfem_interpolation.h
Interpolation of fields from a mesh_fem onto another.
getfem_mesh_slice.h
Define the class getfem::stored_mesh_slice.
gmm::copy
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:976
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:48
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
getfem::interpolate_transformation_neighbor_instance
pinterpolate_transformation interpolate_transformation_neighbor_instance()
Create a new instance of a transformation corresponding to the interpolation on the neighbor element.
Definition: getfem_generic_assembly_interpolation.cc:862
getfem::add_interpolate_transformation_on_deformed_domains
void add_interpolate_transformation_on_deformed_domains(ga_workspace &workspace, const std::string &transname, const mesh &source_mesh, const std::string &source_displacements, const mesh_region &source_region, const mesh &target_mesh, const std::string &target_displacements, const mesh_region &target_region)
Add a transformation to the workspace that creates an identity mapping between two meshes in deformed...
Definition: getfem_interpolation_on_deformed_domains.cc:361
getfem::add_interpolate_transformation_from_expression
void add_interpolate_transformation_from_expression(ga_workspace &workspace, const std::string &transname, const mesh &source_mesh, const mesh &target_mesh, const std::string &expr)
Add a transformation to a workspace workspace or a model md mapping point in mesh source_mesh to mesh...
Definition: getfem_generic_assembly_interpolation.cc:772
getfem::ga_local_projection
void ga_local_projection(const getfem::model &md, const mesh_im &mim, const std::string &expr, const mesh_fem &mf, base_vector &result, const mesh_region &rg=mesh_region::all_convexes())
Make an elementwise L2 projection of an expression with respect to the mesh_fem mf.
Definition: getfem_generic_assembly_interpolation.cc:429
bgeot::stored_point_tab
Point tab storage.
Definition: bgeot_convex_ref.h:49