GetFEM  5.4.2
getfem_config.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2000-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 /*! @mainpage GetFEM reference documentation.
33 
34  <center><img src="http://download.gna.org/getfem/html/homepage/_static/logo_getfem_small.png"></center>
35 
36  @section intro Introduction
37 
38  This documentation has been generated from the GetFEM sources, this is not a tutorial.
39 
40  @section Terminology
41 
42  This is just a short summary of the terms employed in getfem.
43 
44  The @link mesh mesh@endlink is composed of @link convexes
45  convexes@endlink. What we call convexes can be simple line segments,
46  prisms, tetrahedrons, curved triangles, of even something which is
47  not convex (in the geometrical sense). They all have an associated
48  @link refconvexes reference convex@endlink: for segments, this will
49  be the @f$[0,1]@f$ segment, for triangles this will be the canonical
50  triangle @f$(0,0)-(0,1)-(1,0)@f$ etc... All convexes of the mesh
51  are constructed from the reference convex through a @link
52  bgeot_geometric_trans.h geometric transformation@endlink. In
53  simple cases (when the convexes are simplices for example), this
54  transformation will be linear (hence it is easily inverted, which
55  can be a great advantage). In order to define the geometric
56  transformation, one defines <i>geometrical nodes</i> on the
57  reference convex. The geometrical transformation maps these nodes to
58  the <i>mesh nodes</i>.
59 
60  On the mesh, one defines a set a <i>basis functions</i>: the @link pfem FEM@endlink. A FEM
61  should be associated with each convex. The basis functions are also
62  attached to some geometrical points (which can be arbitrarily
63  chosen). These points are similar to the mesh nodes, but <b>they
64  don't have to be the same</b> (this only happens on very simple cases,
65  such as a classical P1 fem on a triangular mesh). The set of all
66  basis functions on the mesh forms the basis of a vector space, on
67  which the PDE will be solved. These basis functions (and their
68  associated geometrical point) are the <b>degrees of freedom
69  (dof)</b>. The FEM is said to be <i>Lagrangian</i> when each of its basis
70  functions is equal to one at its attached geometrical point, and is
71  null at the geometrical points of others basis functions. This is an
72  important property as it is very easy to <i>interpolate</i> an arbitrary
73  function on the finite elements space.
74 
75  The finite elements methods involves evaluation of integrals of
76  these basis functions (or product of basis functions etc...) on
77  convexes (and faces of convexes). In simple cases (polynomial basis
78  functions and linear geometrical transformation), one can evaluate
79  analytically these integrals. In other cases, one has to approximate
80  it, using <i>quadrature formulas</i>. Hence, at each convex is attached
81  an @link getfem_integration.h integration method@endlink along with the FEM. If you have to use an
82  approximate integration method, always choose carefully its
83  order(i.e. highest degree of the polynomials who are exactly
84  integrated with the method) : the degree of the FEM, of the
85  polynomial degree of the geometrical transformation, and the nature
86  of the elementary matrix have to be taken into account. If you are
87  unsure about the appropriate degree, always prefer a high order
88  integration method (which will slow down the assembly) to a low
89  order one which will produce a useless linear-system.
90 
91  The process of construction of a global linear system from integrals
92  of basis functions on each convex is the @link asm assembly@endlink.
93 
94  A mesh, with a set of FEM attached to its convexes is called a @link getfem_mesh_fem.h
95  mesh_fem@endlink object in GetFEM.
96 
97  A mesh, with a set of integration methods attached to its convexes
98  is called a @link getfem_mesh_im.h mesh_im@endlink object in GetFEM.
99 
100  A @c mesh_fem can be used to approximate scalar fields (heat, pression,
101  ..), or vector fields (displacement, electric field, ..). A @c mesh_im
102  will be used to perform numerical integrations on these fields.
103  Although GetFEM supports vector FEMs, intrinsic vector FEMs are
104  essentially used in mixed methods (for instance Raviart-Thomas element).
105  Most of the FEM currently available are scalar. Of course,
106  these scalar FEMs can be used to approximate each component of a
107  vector field. This is done by setting the Qdim of the mesh_fem to
108  the dimension of the vector field (i.e. Qdim=1 for scalar field,
109  Qdim=2 for 2D vector field etc...).
110 
111  When solving a PDE, one often has to use more than one FEM. The most
112  important one will be of course the one on which is defined the
113  solution of the PDE. But most PDEs involve various coefficients, for
114  example: @f[ \nabla.(\lambda(x)\nabla u) = f(x). @f] Hence one has
115  to define an FEM for the main unknown @f$u@f$, but also for the data
116  @f$\lambda(x)@f$ and @f$f(x)@f$ if they are not constant. In order
117  to interpolate easily these coefficients in their finite element
118  space, one often choose a Lagrangian FEM.
119 
120  The convexes, mesh nodes, and dof are all numbered. We sometimes
121  refer to the number associated to a convex as its convex id or
122  convex number. Mesh node numbers are also called point id or point
123  number. Faces of convexes do not have a global numbering, but only
124  a local number in each convex.
125 
126  While the @c dof are always numbered consecutively, <b>this is not
127  always the case for point ids and convex ids</b>, especially if you
128  have removed points or convexes from the mesh. To ensure that they
129  form a continuous sequence (starting from 0), you have to use the
130  getfem::getfem_mesh::optimize_structure member function.
131 
132  Most of the example programs of getfem now uses @link bricks model bricks@endlink.
133  Connecting these basic blocks one to each other, solvers for many PDE problems can be built.
134 
135  @section Examples
136 
137  - @link laplacian.cc tests/laplacian.cc@endlink: solve the laplace equation.
138  - @link elastostatic.cc tests/elastostatic.cc@endlink: solve a static linear elasticity problem.
139  - @link helmholtz.cc tests/helmholtz.cc@endlink: scalar wave equation, with robin boundary conditions.
140  - @link stokes.cc tests/stokes.cc@endlink: the Stokes equation (incompressible viscuous fluid).
141  - @link nonlinear_elastostatic.cc tests/nonlinear_elastostatic.cc@endlink: large strain elastostatic problem (torsion of a bar).
142  - @link icare.cc contrib/icare/icare.cc@endlink: Navier-Stokes equation (fluid flow around an obstacle).
143 */
144 
145 /**@file getfem_config.h
146  @author Yves Renard <Yves.Renard@insa-lyon.fr>
147  @date November 19, 2000.
148  @brief defines and typedefs for namespace getfem
149 */
150 
151 
152 #ifndef GETFEM_CONFIG_H__
153 #define GETFEM_CONFIG_H__
154 
155 #include "bgeot_config.h"
156 
157 // Parallelisation options
158 
159 // GETFEM_PARA_LEVEL is the parallelisation level of Getfem
160 // 0 - Sequential
161 // 1 - Only the resolution of linear systems are parallelized
162 // 2 - Assembly procedures are also parallelized
163 #ifndef GETFEM_PARA_LEVEL
164 # define GETFEM_PARA_LEVEL 0
165 #endif
166 
167 #define GETFEM_MPI_INIT(argc, argv) {GMM_TRACE1("Running sequential Getfem");}
168 #define GETFEM_MPI_FINALIZE {}
169 
170 #if defined(GETFEM_HAVE_DMUMPS_C_H)
171 # ifndef GMM_USES_MUMPS
172 # define GMM_USES_MUMPS
173 # endif
174 #endif
175 
176 
177 #if GMM_USES_MPI > 0
178 # include <mpi.h>
179 
180 # undef GMM_TRACE_MSG_MPI
181 # define GMM_TRACE_MSG_MPI \
182  int mip_rk__; MPI_Comm_rank(MPI_COMM_WORLD, &mip_rk__); \
183  if (mip_rk__ == 0)
184 
185 # undef GETFEM_MPI_INIT
186 # define GETFEM_MPI_INIT(argc, argv) { \
187  MPI_Init(&argc, &argv); \
188  GMM_TRACE1("Running parallelized Getfem level " << GETFEM_PARA_LEVEL); \
189  }
190 # undef GETFEM_MPI_FINALIZE
191 # define GETFEM_MPI_FINALIZE { MPI_Finalize(); }
192 
193 // GETFEM_PARA_SOLVER is the parallelisation solver used
194 // MUMPS : use direct parallel solver MUMPS
195 // SCHWARZADD : use a Schwarz additive method
196 #define MUMPS_PARA_SOLVER 1
197 #define SCHWARZADD_PARA_SOLVER 2
198 
199 # ifndef GETFEM_PARA_SOLVER
200 # define GETFEM_PARA_SOLVER MUMPS_PARA_SOLVER
201 # endif
202 
203 # if GETFEM_PARA_SOLVER == MUMPS_PARA_SOLVER
204 # ifndef GMM_USES_MUMPS
205 # define GMM_USES_MUMPS
206 # endif
207 # endif
208 
209 #endif
210 
211 
212 #include "bgeot_tensor.h"
213 #include "bgeot_poly.h"
214 #include "getfem_superlu.h"
215 
216 /// GEneric Tool for Finite Element Methods.
217 namespace getfem {
218 
219  using std::endl; using std::cout; using std::cerr;
220  using std::ends; using std::cin;
221  using gmm::vref;
222 
223 #if GETFEM_PARA_LEVEL > 1
224  template <typename T> inline T MPI_SUM_SCALAR(T a) {
225  T b; MPI_Allreduce(&a,&b,1,gmm::mpi_type(a), MPI_SUM, MPI_COMM_WORLD);
226  return b;
227  }
228 
229  template <typename VECT> inline void MPI_SUM_VECTOR(const VECT &VV) {
230  VECT &V = const_cast<VECT &>(VV);
231  typedef typename gmm::linalg_traits<VECT>::value_type T;
232  std::vector<T> W(gmm::vect_size(V));
233  MPI_Allreduce((void *)(&(V[0])), &(W[0]), gmm::vect_size(V),
234  gmm::mpi_type(T()), MPI_SUM, MPI_COMM_WORLD);
235  gmm::copy(W, V);
236  }
237 
238  template <typename VECT> inline void MPI_MAX_VECTOR(const VECT &VV) {
239  VECT &V = const_cast<VECT &>(VV);
240  typedef typename gmm::linalg_traits<VECT>::value_type T;
241  std::vector<T> W(gmm::vect_size(V));
242  MPI_Allreduce((void *)(&(V[0])), &(W[0]), gmm::vect_size(V),
243  gmm::mpi_type(T()), MPI_MAX, MPI_COMM_WORLD);
244  gmm::copy(W, V);
245  }
246 
247  template <typename T> void MPI_BCAST0_SCALAR(T &a) {
248  MPI_Bcast((void *)(&a), 1, gmm::mpi_type(a), 0, MPI_COMM_WORLD);
249  }
250 
251  template <typename VECT> inline void MPI_BCAST0_VECTOR(const VECT &VV) {
252  VECT &V = const_cast<VECT &>(VV);
253  typedef typename gmm::linalg_traits<VECT>::value_type T;
254  MPI_Bcast((void *)(&(V[0])), gmm::vect_size(V), gmm::mpi_type(T()), 0,
255  MPI_COMM_WORLD);
256  }
257 
258  template <typename VECT1, typename VECT2>
259  inline void MPI_SUM_VECTOR(const VECT1 &VV, const VECT2 &WW) {
260  VECT1 &V = const_cast<VECT1 &>(VV);
261  VECT2 &W = const_cast<VECT2 &>(WW);
262  typedef typename gmm::linalg_traits<VECT1>::value_type T;
263  MPI_Allreduce((void *)(&(V[0])), &(W[0]),
264  gmm::vect_size(V), gmm::mpi_type(T()),
265  MPI_SUM, MPI_COMM_WORLD);
266  }
267 
268  inline bool MPI_IS_MASTER(void)
269  { int rk; MPI_Comm_rank(MPI_COMM_WORLD, &rk); return !rk; }
270 
271  template <typename T> inline
272  void MPI_SUM_SPARSE_MATRIX(const gmm::col_matrix< gmm::rsvector<T> > &M) {
273  typedef typename gmm::col_matrix< gmm::rsvector<T> > MAT;
274  MAT &MM = const_cast<MAT &>(M);
275  int rk, nbp;
276  MPI_Status status;
277  MPI_Comm_rank(MPI_COMM_WORLD, &rk);
278  MPI_Comm_size(MPI_COMM_WORLD, &nbp);
279  if (nbp < 2) return;
280  size_type nr = gmm::mat_nrows(MM), nc = gmm::mat_ncols(MM);
281 
282  gmm::dense_matrix<int> all_nnz(nc, nbp);
283  std::vector<int> my_nnz(nc), owner(nc);
284  gmm::rsvector<T> V(nr);
285 
286  // Step 1 : each process fill the corresponding nnz line
287  for (size_type i = 0; i < nc; ++i)
288  my_nnz[i] = gmm::nnz(gmm::mat_col(MM, i));
289 
290  // Step 2 : All gather : each proc has all the information
291  MPI_Allgather((void *)(&(my_nnz[0])), nc, MPI_INT,
292  (void *)(&(all_nnz(0,0))), nc, MPI_INT, MPI_COMM_WORLD);
293 
294  // Step 3 : Scan each column and message send/receive
295  std::vector<int> contributors(nc);
296  for (int i = 0; i < nc; ++i) {
297  if (my_nnz[i]) {
298  int column_master = -1, max_nnz = 0;
299  contributors.resize(0);
300 
301  // Determine who is the master for the column
302  for (int j = nbp-1; j >= 0; --j) {
303  if (all_nnz(i, j)) {
304  if (rk != j) contributors.push_back(j);
305  if (column_master == -1 || all_nnz(i, j) > max_nnz)
306  { column_master = j; max_nnz = all_nnz(i, j); }
307  }
308  }
309 
310  if (column_master == rk) { // receive a column and store
311  for (int j = 0; j < int(contributors.size()); ++j) {
312  typename gmm::rsvector<T>::base_type_ &VV = V;
313  int si = all_nnz(i, contributors[j]);
314  VV.resize(si);
315  MPI_Recv((void *)(&(VV[0])),
316  si*sizeof(gmm::elt_rsvector_<T>),
317  MPI_BYTE, contributors[j], 0,
318  MPI_COMM_WORLD, MPI_STATUS_IGNORE);
319  gmm::add(V, gmm::mat_col(MM, i));
320  }
321  } else { // send the column to the column master
322  typename gmm::rsvector<T>::base_type_ &VV = MM.col(i);
323  MPI_Send((void *)(&(VV[0])),
324  my_nnz[i]*sizeof(gmm::elt_rsvector_<T>),
325  MPI_BYTE, column_master, 0,
326  MPI_COMM_WORLD);
327  MM.col(i).clear();
328  }
329  }
330  }
331 
332  // Step 3 : gather the total nnz
333  for (size_type i = 0; i < nc; ++i) {
334  my_nnz[i] = gmm::nnz(gmm::mat_col(MM, i));
335  owner[i] = (my_nnz[i]) ? rk : 0;
336  }
337  MPI_SUM_VECTOR(my_nnz);
338  MPI_SUM_VECTOR(owner);
339 
340  // Step 4 : distribute each column to all the processes
341  // Would it be more efficient to perform a global MPI_SUM on a compressed
342  // storage ?
343  for (size_type i = 0; i < nc; ++i) {
344  if (my_nnz[i]) {
345  typename gmm::rsvector<T>::base_type_ &VV = MM.col(i);
346  if (owner[i] != rk) VV.resize(my_nnz[i]);
347  MPI_Bcast((void *)(&(VV[0])), my_nnz[i]*sizeof(gmm::elt_rsvector_<T>),
348  MPI_BYTE, owner[i], MPI_COMM_WORLD);
349  }
350  }
351  }
352 
353  template <typename MAT> inline void MPI_SUM_SPARSE_MATRIX(const MAT &M) {
354  typedef typename gmm::linalg_traits<MAT>::value_type T;
355  int nbp; MPI_Comm_size(MPI_COMM_WORLD, &nbp);
356  if (nbp < 2) return;
357  size_type nr = gmm::mat_nrows(M), nc = gmm::mat_ncols(M);
358  gmm::col_matrix< gmm::rsvector<T> > MM(nr, nc);
359  GMM_WARNING2("MPI_SUM_SPARSE_MATRIX: A copy of the matrix is done. "
360  "To avoid it, adapt MPI_SUM_SPARSE_MATRIX to "
361  "your matrix type.");
362  gmm::copy(M, MM);
363  MPI_SUM_SPARSE_MATRIX(MM);
364  gmm::copy(MM, const_cast<MAT &>(M));
365  }
366 #else
367  template <typename T> inline T MPI_SUM_SCALAR(T a) { return a; }
368  template <typename VECT> inline void MPI_SUM_VECTOR(const VECT &) {}
369  template <typename VECT> inline void MPI_MAX_VECTOR(const VECT &) {}
370  template <typename T> void MPI_BCAST0_SCALAR(T &) {}
371  template <typename VECT> inline void MPI_BCAST0_VECTOR(const VECT &) {}
372  template <typename MAT> inline void MPI_SUM_SPARSE_MATRIX(const MAT &) {}
373  template <typename VECT1, typename VECT2>
374  inline void MPI_SUM_VECTOR(const VECT1 &V, const VECT2 &WW)
375  {
376  VECT2 &W = const_cast<VECT2 &>(WW);
377  gmm::copy(V, W);
378  }
379  inline bool MPI_IS_MASTER(void) { return true; }
380 #endif
381 
382  using bgeot::ST_NIL;
383  using bgeot::size_type;
384  using bgeot::dim_type;
385  using bgeot::short_type;
386  using bgeot::short_type;
387  using bgeot::scalar_type;
388  using bgeot::complex_type;
389  using bgeot::long_scalar_type;
390  using bgeot::opt_long_scalar_type;
391 
393  using bgeot::base_vector;
394  using bgeot::base_complex_vector;
395  using bgeot::base_matrix;
396  using bgeot::base_complex_matrix;
397  using bgeot::base_tensor;
398  using bgeot::base_complex_tensor;
399  using bgeot::base_poly;
400  using bgeot::base_node;
401 
402 #if defined(__GNUC__)
403  using std::isnan;
404 #else
405  inline bool isnan(scalar_type x) { return x != x; }
406 #endif
407 
408 } /* end of namespace getfem. */
409 
410 
411 #endif /* GETFEM_CONFIG_H__ */
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
getfem_superlu.h
SuperLU interface for getfem.
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
bgeot_tensor.h
tensor class, used in mat_elem computations.
bgeot_poly.h
Multivariate polynomials.
bgeot::polynomial< opt_long_scalar_type >
bgeot::small_vector< scalar_type >
gmm::rsvector
sparse vector built upon std::vector.
Definition: gmm_def.h:488
bgeot_config.h
defines and typedefs for namespace bgeot