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