GetFEM  5.4.2
getfem_convect.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2009-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_convect.h
33  @author Yves Renard <Yves.Renard@insa-lyon.fr>,
34  @date October 27, 2009.
35  @brief Compute the convection of a quantity with respect to a vector field.
36 */
37 #ifndef GETFEM_CONVECT_H__
38 #define GETFEM_CONVECT_H__
39 
40 #include "getfem_mesh_fem.h"
41 #include "getfem_interpolation.h"
42 #include "gmm/gmm_dense_qr.h"
43 
44 namespace getfem {
45 
46  enum convect_boundary_option { CONVECT_EXTRAPOLATION, CONVECT_UNCHANGED, CONVECT_PERIODICITY };
47 
48  /** Compute the convection of a quantity on a getfem::mesh_fem with respect
49  to a velocity field
50  @param mf the source mesh_fem. Should be of Lagrange type.
51  @param U the source field.
52  @param mf_v the mesh_fem on which the vector field is described
53  @param V contains the vector field described on mf_v.
54  @param nt number of time integration step.
55  @param option concerns the entrant boundary.
56  @param per_min, per_max : the periodicity box for the PERIODICITY option
57  */
58  template<class VECT1, class VECT2>
59  void convect(const mesh_fem &mf, VECT1 &U, const mesh_fem &mf_v,
60  const VECT2 &V, scalar_type dt, size_type nt,
61  convect_boundary_option option = CONVECT_EXTRAPOLATION,
62  const base_node &per_min = base_node(),
63  const base_node &per_max = base_node()) {
64  // Should be robustified on the boundaries -> control of the nodes going
65  // out the mesh.
66  // Should control that the point do not move to fast ( < h/2 for instance).
67  // Could be extended to non-lagragian fem with a projection (moving the
68  // gauss points).
69  // Can take into account a source term by integration on the
70  // characteristics.
71 
72  typedef typename gmm::linalg_traits<VECT1>::value_type T;
73  int extra = (option == CONVECT_EXTRAPOLATION) ? 2 : 0;
74 
75  if (nt == 0) return;
76 
77  GMM_ASSERT1(!(mf.is_reduced()),
78  "This convection algorithm work only on pure Lagrange fems");
79  /* test if mf is really of Lagrange type. */
80  for (dal::bv_visitor cv(mf.convex_index()); !cv.finished();++cv) {
81  pfem pf_t = mf.fem_of_element(cv);
82  GMM_ASSERT1(pf_t->target_dim() == 1 && pf_t->is_lagrange(),
83  "This convection algorithm work only on pure Lagrange fems");
84  }
85 
86  // Get the nodes of mf
87  const mesh &msh(mf.linked_mesh());
88  size_type N = msh.dim();
89  if (option == CONVECT_PERIODICITY)
90  GMM_ASSERT1(per_min.size() == N && per_max.size() == N,
91  "Wrong size of box extremity for PERIODICITY option");
92 
93  getfem::mesh_trans_inv mti(msh, 1E-10);
94  size_type qdim = mf.get_qdim();
95  size_type nbpts = mf.nb_basic_dof() / qdim;
96  std::vector<base_node> nodes(nbpts);
97  for (size_type i = 0; i < nbpts; ++i)
98  nodes[i] = mf.point_of_basic_dof(i * qdim);
99 
100  // Obtain the first interpolation of v (same mesh)
101  size_type qqdimt = (gmm::vect_size(V) / mf_v.nb_dof()) * mf_v.get_qdim();
102  GMM_ASSERT1(qqdimt == N, "The velocity field should be a vector field "
103  "of the same dimension as the mesh");
104  std::vector<T> VI(nbpts*N);
105  getfem::interpolation(mf_v, mf, V, VI);
106 
107  // Convect the nodes with respect to v
108  scalar_type ddt = dt / scalar_type(nt);
109  for (size_type i = 0; i < nt; ++i) {
110  if (i > 0) {
111  mti.clear();
112  mti.add_points(nodes);
113  gmm::clear(VI);
114  dal::bit_vector dof_untouched;
115  interpolation(mf_v, mti, V, VI, extra, &dof_untouched);
116  for (dal::bv_visitor j(dof_untouched); !j.finished();++j)
117  VI[j] = V[j];
118  }
119 
120  for (size_type j = 0; j < nbpts; ++j) {
121  gmm::add(gmm::scaled(gmm::sub_vector(VI, gmm::sub_interval(N*j, N)),
122  -ddt), nodes[j]);
123  if (option == CONVECT_PERIODICITY) {
124  for (size_type k = 0; k < N; ++k)
125  if (per_max[k] > per_min[k]) {
126  while (nodes[j][k] > per_max[k]) nodes[j][k] -= per_max[k];
127  while (nodes[j][k] < per_min[k]) nodes[j][k] += per_max[k];
128  }
129  }
130  }
131  }
132 
133  // 3 final interpolation
134  std::vector<T> UI(nbpts*qdim);
135  mti.clear();
136  mti.add_points(nodes);
137  dal::bit_vector dof_untouched;
138  interpolation(mf, mti, U, UI, extra, &dof_untouched);
139  for (dal::bv_visitor i(dof_untouched); !i.finished();++i)
140  UI[i] = U[i];
141  gmm::copy(UI, U);
142  }
143 
144 
145 }
146 
147 
148 #endif
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
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
getfem::convect
void convect(const mesh_fem &mf, VECT1 &U, const mesh_fem &mf_v, const VECT2 &V, scalar_type dt, size_type nt, convect_boundary_option option=CONVECT_EXTRAPOLATION, const base_node &per_min=base_node(), const base_node &per_max=base_node())
Compute the convection of a quantity on a getfem::mesh_fem with respect to a velocity field.
Definition: getfem_convect.h:59
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
Describe a finite element method linked to a mesh.
Definition: getfem_mesh_fem.h:148
getfem::interpolation
void interpolation(const mesh_fem &mf_source, const mesh_fem &mf_target, const VECTU &U, VECTV &V, int extrapolation=0, double EPS=1E-10, mesh_region rg_source=mesh_region::all_convexes(), mesh_region rg_target=mesh_region::all_convexes())
interpolation/extrapolation of (mf_source, U) on mf_target.
Definition: getfem_interpolation.h:693
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
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::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
getfem_mesh_fem.h
Define the getfem::mesh_fem class.
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::linked_mesh
const mesh & linked_mesh() const
Return a reference to the underlying mesh.
Definition: getfem_mesh_fem.h:279
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
Describe a mesh (collection of convexes (elements) and points).
Definition: getfem_mesh.h:95
getfem_interpolation.h
Interpolation of fields from a mesh_fem onto another.
gmm_dense_qr.h
Dense QR factorization.
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