GetFEM  5.4.2
getfem_mesh_fem.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 1999-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_mesh_fem.h
33  @author Yves Renard <Yves.Renard@insa-lyon.fr>
34  @date December 21, 1999.
35  @brief Define the getfem::mesh_fem class
36 */
37 
38 #ifndef GETFEM_MESH_FEM_H__
39 #define GETFEM_MESH_FEM_H__
40 
41 #include "getfem_mesh.h"
42 #include "getfem_fem.h"
43 
44 
45 namespace getfem {
46 
47  template <class CONT> struct tab_scal_to_vect_iterator {
48 
49  typedef typename CONT::const_iterator ITER;
50  typedef typename std::iterator_traits<ITER>::value_type value_type;
51  typedef typename std::iterator_traits<ITER>::pointer pointer;
52  typedef typename std::iterator_traits<ITER>::reference reference;
53  typedef typename std::iterator_traits<ITER>::difference_type
54  difference_type;
55  typedef typename std::iterator_traits<ITER>::iterator_category
56  iterator_category;
57  typedef size_t size_type;
58  typedef tab_scal_to_vect_iterator<CONT> iterator;
59 
60  ITER it;
61  dim_type N;
62  dim_type ii;
63 
64  iterator &operator ++()
65  { ++ii; if (ii == N) { ii = 0; ++it; } return *this; }
66  iterator &operator --()
67  { if (ii == 0) { ii = N; --it; } --ii; return *this; }
68  iterator operator ++(int) { iterator tmp = *this; ++(*this); return tmp; }
69  iterator operator --(int) { iterator tmp = *this; --(*this); return tmp; }
70 
71  iterator &operator +=(difference_type i)
72  { it += (i+ii)/N; ii = dim_type((ii + i) % N); return *this; }
73  iterator &operator -=(difference_type i)
74  { it -= (i+N-ii-1)/N; ii = (ii - i + N * i) % N; return *this; }
75  iterator operator +(difference_type i) const
76  { iterator itt = *this; return (itt += i); }
77  iterator operator -(difference_type i) const
78  { iterator itt = *this; return (itt -= i); }
79  difference_type operator -(const iterator &i) const
80  { return (it - i.it) * N + ii - i.ii; }
81 
82  value_type operator *() const { return (*it) + ii; }
83  value_type operator [](int i) { return *(this + i); }
84 
85  bool operator ==(const iterator &i) const
86  { return (it == i.it) && (ii == i.ii); }
87  bool operator !=(const iterator &i) const { return !(i == *this); }
88  bool operator < (const iterator &i) const
89  { return (it < i.it) && (ii < i.ii); }
90 
91  tab_scal_to_vect_iterator() {}
92  tab_scal_to_vect_iterator(const ITER &iter, dim_type n, dim_type i)
93  : it(iter), N(n), ii(i) { }
94 
95  };
96 
97  /** @internal @brief structure for iteration over the dofs when Qdim
98  != 1 and target_dim == 1
99  */
100  template <class CONT> class tab_scal_to_vect {
101  public :
102  typedef typename CONT::const_iterator ITER;
103  typedef typename std::iterator_traits<ITER>::value_type value_type;
104  typedef typename std::iterator_traits<ITER>::pointer pointer;
105  typedef typename std::iterator_traits<ITER>::pointer const_pointer;
106  typedef typename std::iterator_traits<ITER>::reference reference;
107  typedef typename std::iterator_traits<ITER>::reference const_reference;
108  typedef typename std::iterator_traits<ITER>::difference_type
109  difference_type;
110  typedef size_t size_type;
111  typedef tab_scal_to_vect_iterator<CONT> iterator;
112  typedef iterator const_iterator;
113  typedef std::reverse_iterator<const_iterator> const_reverse_iterator;
114  typedef std::reverse_iterator<iterator> reverse_iterator;
115 
116 
117  protected :
118  ITER it, ite;
119  dim_type N;
120 
121  public :
122 
123  bool empty() const { return it == ite; }
124  size_type size() const { return (ite - it) * N; }
125 
126  const_iterator begin() const { return iterator(it, N, 0); }
127  const_iterator end() const { return iterator(ite, N, 0); }
128  const_reverse_iterator rbegin() const
129  { return const_reverse_iterator(end()); }
130  const_reverse_iterator rend() const
131  { return const_reverse_iterator(begin()); }
132 
133  value_type front() const { return *begin(); }
134  value_type back() const { return *(--(end())); }
135 
136  tab_scal_to_vect() : N(0) {}
137  tab_scal_to_vect(const CONT &cc, dim_type n)
138  : it(cc.begin()), ite(cc.end()), N(n) {}
139 
140  value_type operator [](size_type ii) const { return *(begin() + ii);}
141  };
142 
143  /** Describe a finite element method linked to a mesh.
144  *
145  * @see mesh
146  * @see mesh_im
147  */
148  class mesh_fem : public context_dependencies, virtual public dal::static_stored_object {
149  protected :
150  typedef gmm::csc_matrix<scalar_type> REDUCTION_MATRIX;
151  typedef gmm::csr_matrix<scalar_type> EXTENSION_MATRIX;
152 
153  void copy_from(const mesh_fem &mf); /* Remember to change copy_from if
154  adding components to mesh_fem */
155 
156  std::vector<pfem> f_elems;
157  dal::bit_vector fe_convex;
158  const mesh *linked_mesh_;
159  REDUCTION_MATRIX R_;
160  EXTENSION_MATRIX E_;
161  mutable bgeot::mesh_structure dof_structure;
162  mutable bool dof_enumeration_made;
163  mutable bool is_uniform_, is_uniformly_vectorized_;
164  mutable size_type nb_total_dof;
165  pfem auto_add_elt_pf; /* fem for automatic addition */
166  /* of element option. (0 = no automatic addition) */
167  dim_type auto_add_elt_K; /* Degree of the fem for automatic addition */
168  /* of element option. (-1 = no automatic addition) */
169  bool auto_add_elt_disc, auto_add_elt_complete;
170  scalar_type auto_add_elt_alpha;
171  dim_type Qdim; /* this is the "total" target_dim. */
172  bgeot::multi_index mi; /* multi dimensions for matrix/tensor field. */
173  // dim_type QdimM, QdimN; /* for matrix field with QdimM lines and QdimN */
174  // /* columnsQdimM * QdimN = Qdim. */
175  std::vector<size_type> dof_partition;
176  mutable gmm::uint64_type v_num_update, v_num;
177  bool use_reduction; /* A reduction matrix is applied or not. */
178 
179  public :
180  typedef base_node point_type;
181  typedef tab_scal_to_vect<mesh::ind_cv_ct> ind_dof_ct;
182  typedef tab_scal_to_vect<mesh::ind_pt_face_ct> ind_dof_face_ct;
183 
184  void update_from_context() const;
185 
186  gmm::uint64_type version_number() const
187  { context_check(); return v_num; }
188 
189  /** Get the set of convexes where a finite element has been assigned.
190  */
191  inline const dal::bit_vector &convex_index() const
192  { context_check(); return fe_convex; }
193 
194  /// Return true if a reduction matrix is applied to the dofs.
195  bool is_reduced() const { return use_reduction; }
196 
197  /// Return the reduction matrix applied to the dofs.
198  const REDUCTION_MATRIX &reduction_matrix() const { return R_; }
199 
200  /// Return the extension matrix corresponding to reduction applied (RE=I).
201  const EXTENSION_MATRIX &extension_matrix() const { return E_; }
202 
203  /** Allows to set the reduction and the extension matrices.
204  * Should satify (RR*EE=I). */
205  template <typename MATR, typename MATE>
206  void set_reduction_matrices(const MATR &RR, const MATE &EE) {
207  context_check();
208  GMM_ASSERT1(gmm::mat_ncols(RR) == nb_basic_dof() &&
209  gmm::mat_nrows(EE) == nb_basic_dof() &&
210  gmm::mat_nrows(RR) == gmm::mat_ncols(EE),
211  "Wrong dimension of reduction and/or extension matrices");
212  R_ = REDUCTION_MATRIX(gmm::mat_nrows(RR), gmm::mat_ncols(RR));
213  E_ = EXTENSION_MATRIX(gmm::mat_nrows(EE), gmm::mat_ncols(EE));
214  gmm::copy(RR, R_);
215  gmm::copy(EE, E_);
216  use_reduction = true;
217  touch(); v_num = act_counter();
218  }
219 
220  /** Allows to set the reduction and the extension matrices in order
221  * to keep only a certain number of dof. */
222  void reduce_to_basic_dof(const dal::bit_vector &kept_basic_dof);
223  void reduce_to_basic_dof(const std::set<size_type> &kept_basic_dof);
224 
225  /// Validate or invalidate the reduction (keeping the reduction matrices).
226  void set_reduction(bool r) {
227  if (r != use_reduction) {
228  use_reduction = r;
229  if (use_reduction) {
230  context_check();
231  GMM_ASSERT1(gmm::mat_ncols(R_) == nb_basic_dof() &&
232  gmm::mat_nrows(E_) == nb_basic_dof() &&
233  gmm::mat_nrows(R_) == gmm::mat_ncols(E_),
234  "Wrong dimension of reduction and/or extension matrices");
235  }
236  touch(); v_num = act_counter();
237  }
238  }
239 
240  template<typename VECT1, typename VECT2>
241  void reduce_vector(const VECT1 &V1, const VECT2 &V2) const {
242  if (is_reduced()) {
243  size_type qqdim = gmm::vect_size(V1) / nb_basic_dof();
244  if (qqdim == 1)
245  gmm::mult(reduction_matrix(), V1, const_cast<VECT2 &>(V2));
246  else
247  for (size_type k = 0; k < qqdim; ++k)
248  gmm::mult(reduction_matrix(),
249  gmm::sub_vector(V1, gmm::sub_slice(k, nb_basic_dof(),
250  qqdim)),
251  gmm::sub_vector(const_cast<VECT2 &>(V2),
252  gmm::sub_slice(k, nb_dof(),
253  qqdim)));
254  }
255  else gmm::copy(V1, const_cast<VECT2 &>(V2));
256  }
257 
258  template<typename VECT1, typename VECT2>
259  void extend_vector(const VECT1 &V1, const VECT2 &V2) const {
260  size_type nbd = nb_dof();
261  if (is_reduced() && nbd) {
262  size_type qqdim = gmm::vect_size(V1) / nbd;
263  if (qqdim == 1)
264  gmm::mult(extension_matrix(), V1, const_cast<VECT2 &>(V2));
265  else
266  for (size_type k = 0; k < qqdim; ++k)
267  gmm::mult(extension_matrix(),
268  gmm::sub_vector(V1, gmm::sub_slice(k, nb_dof(),
269  qqdim)),
270  gmm::sub_vector(const_cast<VECT2 &>(V2),
271  gmm::sub_slice(k, nb_basic_dof(),
272  qqdim)));
273  }
274  else gmm::copy(V1, const_cast<VECT2 &>(V2));
275  }
276 
277 
278  /// Return a reference to the underlying mesh.
279  const mesh &linked_mesh() const { return *linked_mesh_; }
280 
281  virtual bool is_uniform() const;
282  virtual bool is_uniformly_vectorized() const;
283 
284  /** Set the degree of the fem for automatic addition
285  * of element option. K=-1 disables the automatic addition.
286  */
287  void set_auto_add(dim_type K, bool disc = false,
288  scalar_type alpha = scalar_type(0),
289  bool complete=false) {
290  auto_add_elt_K = K;
291  auto_add_elt_disc = disc;
292  auto_add_elt_alpha = alpha;
293  auto_add_elt_complete = complete;
294  auto_add_elt_pf = 0;
295  }
296 
297  /** Set the fem for automatic addition
298  * of element option. pf=0 disables the automatic addition.
299  */
300  void set_auto_add(pfem pf) {
301  auto_add_elt_pf = pf;
302  auto_add_elt_K = dim_type(-1);
303  auto_add_elt_disc = false;
304  auto_add_elt_alpha = scalar_type(0);
305  auto_add_elt_complete = false;
306  }
307 
308  /** Return the Q dimension. A mesh_fem used for scalar fields has
309  Q=1, for vector fields, Q is typically equal to
310  linked_mesh().dim().
311  */
312  virtual dim_type get_qdim() const { return Qdim; }
313  virtual const bgeot::multi_index &get_qdims() const { return mi; }
314 
315  /** Change the Q dimension */
316  virtual void set_qdim(dim_type q) {
317  if (q != Qdim || mi.size() != 1) {
318  mi.resize(1); mi[0] = q; Qdim = q;
319  dof_enumeration_made = false; touch(); v_num = act_counter();
320  }
321  }
322 
323  /** Set the dimension for a matrix field. */
324  virtual void set_qdim(dim_type M, dim_type N) {
325  if (mi.size() != 2 || mi[0] != M || mi[1] != N) {
326  mi.resize(2); mi[0] = M; mi[1] = N; Qdim = dim_type(M*N);
327  dof_enumeration_made = false; touch(); v_num = act_counter();
328  }
329  }
330 
331  /** Set the dimension for a fourth order tensor field. */
332  virtual void set_qdim(dim_type M, dim_type N, dim_type O, dim_type P) {
333  if (mi.size() != 4 || mi[0] != M || mi[1] != N || mi[2] != O
334  || mi[3] != P) {
335  mi.resize(4); mi[0] = M; mi[1] = N; mi[2] = O; mi[3] = P;
336  Qdim = dim_type(M*N*O*P);
337  dof_enumeration_made = false; touch(); v_num = act_counter();
338  }
339  }
340 
341  /** Set the dimension for an arbitrary order tensor field. */
342  virtual void set_qdim(const bgeot::multi_index &mii) {
343  GMM_ASSERT1(mii.size() < 7,
344  "Tensor field are taken into account up to order 6.");
345  GMM_ASSERT1(mi.size(), "Wrong sizes");
346  if (!(mi.is_equal(mii))) {
347  mi = mii;
348  Qdim = dim_type(1);
349  for (size_type i = 0; i < mi.size(); ++i)
350  Qdim = dim_type(Qdim*mi[i]);
351  GMM_ASSERT1(Qdim, "Wrong sizes");
352  dof_enumeration_made = false; touch(); v_num = act_counter();
353  }
354  }
355 
356 
357  void set_qdim_mn(dim_type M, dim_type N) IS_DEPRECATED
358  { set_qdim(M,N); }
359 
360  /** for matrix fields, return the number of rows. */
361  dim_type get_qdim_m() const IS_DEPRECATED { return dim_type(mi[0]); }
362  /** for matrix fields, return the number of columns. */
363  dim_type get_qdim_n() const IS_DEPRECATED { return dim_type(mi[1]); }
364 
365 
366  /** Set the finite element method of a convex.
367  @param cv the convex number.
368  @param pf the finite element.
369  */
370  void set_finite_element(size_type cv, pfem pf);
371  /** Set the finite element on a set of convexes.
372  @param cvs the set of convex indexes, as a dal::bit_vector.
373  @param pf the finite element, typically obtained with
374  @code getfem::fem_descriptor("FEM_SOMETHING(..)")
375  @endcode
376  */
377  void set_finite_element(const dal::bit_vector &cvs, pfem pf);
378  /** shortcut for set_finite_element(linked_mesh().convex_index(),pf);
379  and set_auto_add(pf). */
380  void set_finite_element(pfem pf);
381  /** Set a classical (i.e. lagrange polynomial) finite element on
382  a convex.
383  @param cv is the convex number.
384  @param fem_degree the polynomial degree of the finite element.
385  @param complete is a flag for excluding incomplete element
386  irrespective of the element geometric transformation.
387  */
388  void set_classical_finite_element(size_type cv, dim_type fem_degree,
389  bool complete=false);
390  /** Set a classical (i.e. lagrange polynomial) finite element on
391  a set of convexes.
392  @param cvs the set of convexes, as a dal::bit_vector.
393  @param fem_degree the polynomial degree of the finite element.
394  */
395  void set_classical_finite_element(const dal::bit_vector &cvs,
396  dim_type fem_degree,
397  bool complete=false);
398  /** Similar to set_classical_finite_element, but uses
399  discontinuous lagrange elements.
400 
401  @param cv the convex index.
402  @param fem_degree the polynomial degree of the finite element.
403  @param alpha the node inset, 0 <= alpha < 1, where 0 implies
404  usual dof nodes, greater values move the nodes toward the
405  center of gravity, and 1 means that all degrees of freedom
406  collapse on the center of gravity.
407  */
409  dim_type fem_degree,
410  scalar_type alpha=0,
411  bool complete=false);
412  /** Similar to set_classical_finite_element, but uses
413  discontinuous lagrange elements.
414 
415  @param cvs the set of convexes, as a dal::bit_vector.
416  @param fem_degree the polynomial degree of the finite element.
417  @param alpha the node inset, 0 <= alpha < 1, where 0 implies
418  usual dof nodes, greater values move the nodes toward the
419  center of gravity, and 1 means that all degrees of freedom
420  collapse on the center of gravity.
421  */
422  void set_classical_discontinuous_finite_element(const dal::bit_vector &cvs,
423  dim_type fem_degree,
424  scalar_type alpha=0,
425  bool complete=false);
426  /** Shortcut for
427  * set_classical_finite_element(linked_mesh().convex_index(),...)
428  */
429  void set_classical_finite_element(dim_type fem_degree,
430  bool complete=false);
431  /** Shortcut for
432  * set_classical_discontinuous_finite_element(linked_mesh().convex_index()
433  * ,...)
434  */
435  void set_classical_discontinuous_finite_element(dim_type fem_degree,
436  scalar_type alpha=0,
437  bool complete=false);
438  /** Return the basic fem associated with an element (if no fem is
439  * associated, the function will crash! use the convex_index() of
440  * the mesh_fem to check that a fem is associated to a given
441  * convex). This fem does not take into account the optional
442  * vectorization due to qdim nor the optional reduction.
443  */
444  virtual pfem fem_of_element(size_type cv) const
445  { return ((cv < f_elems.size()) ? f_elems[cv] : 0); }
446  /** Give an array of the dof numbers a of convex.
447  * @param cv the convex number.
448  * @return a pseudo-container of the dof number.
449  */
450  virtual ind_dof_ct ind_basic_dof_of_element(size_type cv) const {
451  context_check(); if (!dof_enumeration_made) enumerate_dof();
452  return ind_dof_ct(dof_structure.ind_points_of_convex(cv),
453  dim_type(Qdim/fem_of_element(cv)->target_dim()));
454  }
455  ind_dof_ct ind_dof_of_element(size_type cv) const IS_DEPRECATED
456  { return ind_basic_dof_of_element(cv); }
457  virtual const std::vector<size_type> &
458  ind_scalar_basic_dof_of_element(size_type cv) const
459  { return dof_structure.ind_points_of_convex(cv); }
460  /** Give an array of the dof numbers lying of a convex face (all
461  degrees of freedom whose associated base function is non-zero
462  on the convex face).
463  @param cv the convex number.
464  @param f the face number.
465  @return a pseudo-container of the dof number.
466  */
467  virtual ind_dof_face_ct
469  context_check(); if (!dof_enumeration_made) enumerate_dof();
470  return ind_dof_face_ct
471  (dof_structure.ind_points_of_face_of_convex(cv, f),
472  dim_type(Qdim/fem_of_element(cv)->target_dim()));
473  }
474  ind_dof_face_ct
475  ind_dof_of_face_of_element(size_type cv,short_type f) const IS_DEPRECATED
476  { return ind_basic_dof_of_face_of_element(cv, f); }
477  /** Return the number of dof lying on the given convex face.
478  @param cv the convex number.
479  @param f the face number.
480  */
482  short_type f) const {
483  context_check(); if (!dof_enumeration_made) enumerate_dof();
484  pfem pf = f_elems[cv];
485  return dof_structure.structure_of_convex(cv)->nb_points_of_face(f)
486  * Qdim / pf->target_dim();
487  }
488  size_type nb_dof_of_face_of_element(size_type cv,
489  short_type f) const IS_DEPRECATED
490  { return nb_basic_dof_of_face_of_element(cv, f); }
491  /** Return the number of degrees of freedom attached to a given convex.
492  @param cv the convex number.
493  */
495  context_check(); if (!dof_enumeration_made) enumerate_dof();
496  pfem pf = f_elems[cv]; return pf->nb_dof(cv) * Qdim / pf->target_dim();
497  }
498  size_type nb_dof_of_element(size_type cv) const IS_DEPRECATED
499  { return nb_basic_dof_of_element(cv); }
500 
501  /* Return the geometrical location of a degree of freedom in the
502  reference convex.
503  @param cv the convex number.
504  @param i the local dof number.
505  const base_node &reference_point_of_dof(size_type cv,size_type i) const {
506  pfem pf = f_elems[cv];
507  return pf->node_of_dof(cv, i * pf->target_dim() / Qdim);
508  }
509  */
510  /** Return the geometrical location of a degree of freedom.
511  @param cv the convex number.
512  @param i the local dof number.
513  */
514  virtual base_node point_of_basic_dof(size_type cv, size_type i) const;
515  base_node point_of_dof(size_type cv, size_type i) const IS_DEPRECATED
516  { return point_of_basic_dof(cv, i); }
517  /** Return the geometrical location of a degree of freedom.
518  @param d the global dof number.
519  */
520  virtual base_node point_of_basic_dof(size_type d) const;
521  base_node point_of_dof(size_type d) const IS_DEPRECATED
522  { return point_of_basic_dof(d); }
523  /** Return the dof component number (0<= x <Qdim) */
524  virtual dim_type basic_dof_qdim(size_type d) const;
525  dim_type dof_qdim(size_type d) const IS_DEPRECATED
526  { return basic_dof_qdim(d); }
527  /** Shortcut for convex_to_dof(d)[0]
528  @param d the global dof number.
529  */
531  size_type first_convex_of_dof(size_type d) const IS_DEPRECATED
532  { return first_convex_of_basic_dof(d); }
533  /** Return the list of convexes attached to the specified dof
534  @param d the global dof number.
535  @return an array of convex numbers.
536  */
537  virtual const mesh::ind_cv_ct &convex_to_basic_dof(size_type d) const;
538  const mesh::ind_cv_ct &convex_to_dof(size_type d) const IS_DEPRECATED
539  { return convex_to_basic_dof(d); }
540 
541  /** Give an array that contains the global dof indices corresponding
542  * to the mesh_fem dofs or size_type(-1) if a dof is not global.
543  * @param ind the returned global dof indices array.
544  */
545  virtual void get_global_dof_index(std::vector<size_type> &ind) const;
546  /** Renumber the degrees of freedom. You should not have
547  * to call this function, as it is done automatically */
548  virtual void enumerate_dof() const;
549 
550 #if GETFEM_PARA_LEVEL > 1
551  void enumerate_dof_para()const;
552 #endif
553 
554  /** Return the total number of basic degrees of freedom (before the
555  * optional reduction). */
556  virtual size_type nb_basic_dof() const {
557  context_check();
558  if (!dof_enumeration_made) enumerate_dof();
559  return nb_total_dof;
560  }
561  /// Return the total number of degrees of freedom.
562  virtual size_type nb_dof() const {
563  context_check();
564  if (!dof_enumeration_made) enumerate_dof();
565  return use_reduction ? gmm::mat_nrows(R_) : nb_total_dof;
566  }
567  /** Get a list of basic dof lying on a given mesh_region.
568  @param b the mesh_region.
569  @return the list in a dal::bit_vector.
570  */
571  virtual dal::bit_vector basic_dof_on_region(const mesh_region &b) const;
572  /** Get a list of dof lying on a given mesh_region. For a reduced mesh_fem
573  a dof is lying on a region if its potential corresponding shape
574  function is nonzero on this region. The extension matrix is used
575  to make the correspondence between basic and reduced dofs.
576  @param b the mesh_region.
577  @return the list in a dal::bit_vector.
578  */
579  dal::bit_vector dof_on_region(const mesh_region &b) const;
580  dal::bit_vector dof_on_set(const mesh_region &b) const IS_DEPRECATED
581  { return dof_on_region(b); }
582 
583  void set_dof_partition(size_type cv, unsigned partition_num) {
584  if (dof_partition.empty() && partition_num == 0) return;
585  if (dof_partition.size() < linked_mesh().nb_allocated_convex())
586  dof_partition.resize(linked_mesh().nb_allocated_convex());
587  if (dof_partition.at(cv) != partition_num) {
588  dof_partition[cv] = partition_num;
589  dof_enumeration_made = false;
590  }
591  }
592  unsigned get_dof_partition(size_type cv) const {
593  return (cv < dof_partition.size() ? unsigned(dof_partition[cv]) : 0);
594  }
595  void clear_dof_partition() { dof_partition.clear(); }
596 
597  size_type memsize() const {
598  return dof_structure.memsize() +
599  sizeof(mesh_fem) - sizeof(bgeot::mesh_structure) +
600  f_elems.size() * sizeof(pfem) + fe_convex.memsize();
601  }
602  void init_with_mesh(const mesh &me, dim_type Q = 1);
603  /** Build a new mesh_fem. A mesh object must be supplied.
604  @param me the linked mesh.
605  @param Q the Q dimension (see mesh_fem::get_qdim).
606  */
607  explicit mesh_fem(const mesh &me, dim_type Q = 1);
608  mesh_fem();
609  mesh_fem(const mesh_fem &mf);
610  mesh_fem &operator=(const mesh_fem &mf);
611 
612  virtual ~mesh_fem();
613  virtual void clear();
614  /** Read the mesh_fem from a stream.
615  @param ist the stream.
616  */
617  virtual void read_from_file(std::istream &ist);
618  /** Read the mesh_fem from a file.
619  @param name the file name. */
620  void read_from_file(const std::string &name);
621  /* internal usage. */
622  void write_basic_to_file(std::ostream &ost) const;
623  /* internal usage. */
624  void write_reduction_matrices_to_file(std::ostream &ost) const;
625  /** Write the mesh_fem to a stream. */
626  virtual void write_to_file(std::ostream &ost) const;
627  /** Write the mesh_fem to a file.
628 
629  @param name the file name
630 
631  @param with_mesh if set, then the linked_mesh() will also be
632  saved to the file.
633  */
634  void write_to_file(const std::string &name, bool with_mesh=false) const;
635  };
636 
637  /** Gives the descriptor of a classical finite element method of degree K
638  on mesh.
639 
640  The mesh_fem won't be destroyed until its linked_mesh is
641  destroyed. All the mesh_fem built by this function are stored
642  in a cache, which means that calling this function twice with
643  the same arguments will return the same mesh_fem object. A
644  consequence is that you should NEVER modify this mesh_fem!
645  */
646  const mesh_fem &classical_mesh_fem(const mesh &mesh, dim_type degree,
647  dim_type qdim=1, bool complete=false);
648 
649  /** Dummy mesh_fem for default parameter of functions. */
650  const mesh_fem &dummy_mesh_fem();
651 
652 
653  /** Given a mesh_fem @param mf and a vector @param vec of size equal to
654  * mf.nb_basic_dof(), the output vector @param coeff will contain the
655  * values of @param vec corresponding to the basic dofs of element
656  * @param cv. The size of @param coeff is adjusted if necessary.
657  */
658  template <typename VEC1, typename VEC2>
660  const VEC1 &vec,
661  size_type cv, VEC2 &coeff,
662  size_type qmult1 = size_type(-1),
663  size_type qmult2 = size_type(-1)) {
664  if (qmult1 == size_type(-1)) {
665  size_type nbdof = mf.nb_basic_dof();
666  qmult1 = gmm::vect_size(vec) / nbdof;
667  GMM_ASSERT1(gmm::vect_size(vec) == qmult1 * nbdof, "Bad dof vector size");
668  }
669  if (qmult2 == size_type(-1)) {
670  qmult2 = mf.get_qdim();
671  if (qmult2 > 1) qmult2 /= mf.fem_of_element(cv)->target_dim();
672  }
673  size_type qmultot = qmult1*qmult2;
674  auto &ct = mf.ind_scalar_basic_dof_of_element(cv);
675  gmm::resize(coeff, ct.size()*qmultot);
676 
677  auto it = ct.begin();
678  auto itc = coeff.begin();
679  if (qmultot == 1) {
680  for (; it != ct.end(); ++it) *itc++ = vec[*it];
681  } else {
682  for (; it != ct.end(); ++it) {
683  auto itv = vec.begin()+(*it)*qmult1;
684  for (size_type m = 0; m < qmultot; ++m) *itc++ = *itv++;
685  }
686  }
687  }
688 
689  void vectorize_base_tensor(const base_tensor &t, base_matrix &vt,
690  size_type ndof, size_type qdim, size_type N);
691 
692  void vectorize_grad_base_tensor(const base_tensor &t, base_tensor &vt,
693  size_type ndof, size_type qdim, size_type N);
694 
695 
696 } /* end of namespace getfem. */
697 
698 
699 #endif /* GETFEM_MESH_FEM_H__ */
getfem::mesh_fem::enumerate_dof
virtual void enumerate_dof() const
Renumber the degrees of freedom.
Definition: getfem_mesh_fem.cc:321
getfem::mesh_fem::get_qdim
virtual dim_type get_qdim() const
Return the Q dimension.
Definition: getfem_mesh_fem.h:312
getfem::mesh_fem::is_reduced
bool is_reduced() const
Return true if a reduction matrix is applied to the dofs.
Definition: getfem_mesh_fem.h:195
getfem::mesh_fem::set_qdim
virtual void set_qdim(const bgeot::multi_index &mii)
Set the dimension for an arbitrary order tensor field.
Definition: getfem_mesh_fem.h:342
getfem::mesh_fem::convex_to_basic_dof
virtual const mesh::ind_cv_ct & convex_to_basic_dof(size_type d) const
Return the list of convexes attached to the specified dof.
Definition: getfem_mesh_fem.cc:267
getfem::mesh_fem::first_convex_of_basic_dof
virtual size_type first_convex_of_basic_dof(size_type d) const
Shortcut for convex_to_dof(d)[0].
Definition: getfem_mesh_fem.cc:258
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
getfem::mesh_fem::reduce_to_basic_dof
void reduce_to_basic_dof(const dal::bit_vector &kept_basic_dof)
Allows to set the reduction and the extension matrices in order to keep only a certain number of dof.
Definition: getfem_mesh_fem.cc:450
getfem::mesh_fem::set_finite_element
void set_finite_element(size_type cv, pfem pf)
Set the finite element method of a convex.
Definition: getfem_mesh_fem.cc:127
getfem::mesh_fem::basic_dof_qdim
virtual dim_type basic_dof_qdim(size_type d) const
Return the dof component number (0<= x <Qdim)
Definition: getfem_mesh_fem.cc:245
getfem::mesh_fem::dof_on_region
dal::bit_vector dof_on_region(const mesh_region &b) const
Get a list of dof lying on a given mesh_region.
Definition: getfem_mesh_fem.cc:114
getfem::mesh_fem::fem_of_element
virtual pfem fem_of_element(size_type cv) const
Return the basic fem associated with an element (if no fem is associated, the function will crash!...
Definition: getfem_mesh_fem.h:444
getfem::mesh_fem::read_from_file
virtual void read_from_file(std::istream &ist)
Read the mesh_fem from a stream.
Definition: getfem_mesh_fem.cc:540
getfem::mesh_fem::update_from_context
void update_from_context() const
this function has to be defined and should update the object when the context is modified.
Definition: getfem_mesh_fem.cc:30
getfem::mesh_fem::set_auto_add
void set_auto_add(dim_type K, bool disc=false, scalar_type alpha=scalar_type(0), bool complete=false)
Set the degree of the fem for automatic addition of element option.
Definition: getfem_mesh_fem.h:287
bgeot::mesh_structure::ind_points_of_convex
const ind_set & ind_points_of_convex(size_type ic) const
Return a container to the list of points attached to convex ic.
Definition: bgeot_mesh_structure.h:107
getfem::mesh_fem::set_qdim
virtual void set_qdim(dim_type M, dim_type N)
Set the dimension for a matrix field.
Definition: getfem_mesh_fem.h:324
getfem::mesh_fem::get_qdim_n
dim_type get_qdim_n() const IS_DEPRECATED
for matrix fields, return the number of columns.
Definition: getfem_mesh_fem.h:363
getfem_mesh.h
Define a getfem::getfem_mesh object.
getfem::mesh_fem
Describe a finite element method linked to a mesh.
Definition: getfem_mesh_fem.h:148
getfem::mesh_fem::extension_matrix
const EXTENSION_MATRIX & extension_matrix() const
Return the extension matrix corresponding to reduction applied (RE=I).
Definition: getfem_mesh_fem.h:201
getfem::dummy_mesh_fem
const mesh_fem & dummy_mesh_fem()
Dummy mesh_fem for default parameter of functions.
Definition: getfem_mesh_fem.cc:875
getfem::mesh_fem::ind_basic_dof_of_face_of_element
virtual ind_dof_face_ct ind_basic_dof_of_face_of_element(size_type cv, short_type f) const
Give an array of the dof numbers lying of a convex face (all degrees of freedom whose associated base...
Definition: getfem_mesh_fem.h:468
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_fem::basic_dof_on_region
virtual dal::bit_vector basic_dof_on_region(const mesh_region &b) const
Get a list of basic dof lying on a given mesh_region.
Definition: getfem_mesh_fem.cc:77
getfem::mesh_fem::convex_index
const dal::bit_vector & convex_index() const
Get the set of convexes where a finite element has been assigned.
Definition: getfem_mesh_fem.h:191
getfem::mesh_fem::set_qdim
virtual void set_qdim(dim_type M, dim_type N, dim_type O, dim_type P)
Set the dimension for a fourth order tensor field.
Definition: getfem_mesh_fem.h:332
getfem::mesh_fem::nb_basic_dof_of_element
virtual size_type nb_basic_dof_of_element(size_type cv) const
Return the number of degrees of freedom attached to a given convex.
Definition: getfem_mesh_fem.h:494
getfem::mesh_fem::nb_basic_dof
virtual size_type nb_basic_dof() const
Return the total number of basic degrees of freedom (before the optional reduction).
Definition: getfem_mesh_fem.h:556
bgeot::operator+
rational_fraction< T > operator+(const polynomial< T > &P, const rational_fraction< T > &Q)
Add Q to P.
Definition: bgeot_poly.h:749
getfem::context_dependencies
Deal with interdependencies of objects.
Definition: getfem_context.h:81
getfem::mesh_fem::write_to_file
virtual void write_to_file(std::ostream &ost) const
Write the mesh_fem to a stream.
Definition: getfem_mesh_fem.cc:777
bgeot::mesh_structure::ind_points_of_face_of_convex
ind_pt_face_ct ind_points_of_face_of_convex(size_type ic, short_type f) const
Return a container of the (global) point number for face f or convex ic.
Definition: bgeot_mesh_structure.cc:173
getfem::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:244
getfem::mesh_fem::get_global_dof_index
virtual void get_global_dof_index(std::vector< size_type > &ind) const
Give an array that contains the global dof indices corresponding to the mesh_fem dofs or size_type(-1...
Definition: getfem_mesh_fem.cc:293
getfem::mesh_fem::set_reduction_matrices
void set_reduction_matrices(const MATR &RR, const MATE &EE)
Allows to set the reduction and the extension matrices.
Definition: getfem_mesh_fem.h:206
getfem::mesh_fem::linked_mesh
const mesh & linked_mesh() const
Return a reference to the underlying mesh.
Definition: getfem_mesh_fem.h:279
getfem::mesh_fem::set_reduction
void set_reduction(bool r)
Validate or invalidate the reduction (keeping the reduction matrices).
Definition: getfem_mesh_fem.h:226
bgeot::mesh_structure
Mesh structure definition.
Definition: bgeot_mesh_structure.h:73
getfem::mesh_fem::get_qdim_m
dim_type get_qdim_m() const IS_DEPRECATED
for matrix fields, return the number of rows.
Definition: getfem_mesh_fem.h:361
getfem::mesh_fem::point_of_basic_dof
virtual base_node point_of_basic_dof(size_type cv, size_type i) const
Return the geometrical location of a degree of freedom.
Definition: getfem_mesh_fem.cc:223
getfem::mesh_fem::mesh_fem
mesh_fem(const mesh &me, dim_type Q=1)
Build a new mesh_fem.
Definition: getfem_mesh_fem.cc:528
dal::static_stored_object
base class for static stored objects
Definition: dal_static_stored_objects.h:206
getfem::slice_vector_on_basic_dof_of_element
void slice_vector_on_basic_dof_of_element(const mesh_fem &mf, const VEC1 &vec, size_type cv, VEC2 &coeff, size_type qmult1=size_type(-1), size_type qmult2=size_type(-1))
Given a mesh_fem.
Definition: getfem_mesh_fem.h:659
getfem::mesh
Describe a mesh (collection of convexes (elements) and points).
Definition: getfem_mesh.h:95
getfem::mesh_fem::ind_basic_dof_of_element
virtual ind_dof_ct ind_basic_dof_of_element(size_type cv) const
Give an array of the dof numbers a of convex.
Definition: getfem_mesh_fem.h:450
getfem::mesh_region
structure used to hold a set of convexes and/or convex faces.
Definition: getfem_mesh_region.h:55
getfem_fem.h
Definition of the finite element methods.
getfem::context_dependencies::context_check
bool context_check() const
return true if update_from_context was called
Definition: getfem_context.h:126
bgeot::operator-
rational_fraction< T > operator-(const polynomial< T > &P, const rational_fraction< T > &Q)
Subtract Q from P.
Definition: bgeot_poly.h:756
getfem::mesh_fem::set_auto_add
void set_auto_add(pfem pf)
Set the fem for automatic addition of element option.
Definition: getfem_mesh_fem.h:300
getfem::mesh_fem::set_classical_finite_element
void set_classical_finite_element(size_type cv, dim_type fem_degree, bool complete=false)
Set a classical (i.e.
Definition: getfem_mesh_fem.cc:174
getfem::mesh_fem::set_classical_discontinuous_finite_element
void set_classical_discontinuous_finite_element(size_type cv, dim_type fem_degree, scalar_type alpha=0, bool complete=false)
Similar to set_classical_finite_element, but uses discontinuous lagrange elements.
Definition: getfem_mesh_fem.cc:200
getfem::mesh_fem::reduction_matrix
const REDUCTION_MATRIX & reduction_matrix() const
Return the reduction matrix applied to the dofs.
Definition: getfem_mesh_fem.h:198
getfem::mesh_fem::nb_dof
virtual size_type nb_dof() const
Return the total number of degrees of freedom.
Definition: getfem_mesh_fem.h:562
getfem::mesh_fem::set_qdim
virtual void set_qdim(dim_type q)
Change the Q dimension.
Definition: getfem_mesh_fem.h:316
getfem::mesh_fem::nb_basic_dof_of_face_of_element
virtual size_type nb_basic_dof_of_face_of_element(size_type cv, short_type f) const
Return the number of dof lying on the given convex face.
Definition: getfem_mesh_fem.h:481
bgeot::mesh_structure::structure_of_convex
pconvex_structure structure_of_convex(size_type ic) const
Return the pconvex_structure of the convex ic.
Definition: bgeot_mesh_structure.h:112
getfem::classical_mesh_fem
const mesh_fem & classical_mesh_fem(const mesh &mesh, dim_type degree, dim_type qdim=1, bool complete=false)
Gives the descriptor of a classical finite element method of degree K on mesh.
Definition: getfem_mesh_fem.cc:862