GetFEM  5.4.2
getfem_interpolation.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2001-2020 Yves Renard, Julien Pommier
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_interpolation.h
33  @author Yves Renard <[email protected]>
34  @author Julien Pommier <[email protected]>
35  @date October 15, 2001.
36  @brief Interpolation of fields from a mesh_fem onto another.
37 */
38 
39 #ifndef GETFEM_INTERPOLATION_H__
40 #define GETFEM_INTERPOLATION_H__
41 
42 #include "getfem_mesh_fem.h"
43 #include "bgeot_torus.h"
44 #include "dal_tree_sorted.h"
45 #include "getfem_im_data.h"
46 #include "getfem_torus.h"
47 
48 namespace getfem {
49 
50  /* ********************************************************************* */
51  /* */
52  /* I. Distribution of a set of points on a mesh. */
53  /* */
54  /* ********************************************************************* */
55 
56  class mesh_trans_inv : public bgeot::geotrans_inv {
57 
58  protected :
59  typedef std::set<size_type>::const_iterator set_iterator;
60  typedef std::map<size_type,size_type>::const_iterator map_iterator;
61 
62  const mesh &msh;
63  std::vector<std::set<size_type> > pts_cvx;
64  std::vector<base_node> ref_coords;
65  std::map<size_type,size_type> ids;
66 
67  public :
68 
69  size_type nb_points_on_convex(size_type i) const
70  { return pts_cvx[i].size(); }
71  void points_on_convex(size_type i, std::vector<size_type> &itab) const;
72  size_type point_on_convex(size_type cv, size_type i) const;
73  const std::vector<base_node> &reference_coords(void) const { return ref_coords; }
74 
75  void add_point_with_id(base_node n, size_type id)
76  { size_type ipt = add_point(n); ids[ipt] = id; }
77  size_type id_of_point(size_type ipt) const;
78  const mesh &linked_mesh(void) const { return msh; }
79 
80  /* extrapolation = 0 : Only the points inside the mesh are distributed.
81  * extrapolation = 1 : Try to extrapolate the exterior points near the
82  * boundary.
83  * extrapolation = 2 : Extrapolate all the exterior points. Could be
84  * expensive.
85  *
86  * if rg_source is provided only the corresponding part of the mesh is
87  * taken into account and extrapolation is done with respect to the
88  * boundary of the specified region. rg_source must contain only convexes.
89  */
90  void distribute(int extrapolation = 0,
91  mesh_region rg_source=mesh_region::all_convexes());
92  mesh_trans_inv(const mesh &m, double EPS_ = 1E-12)
93  : bgeot::geotrans_inv(EPS_), msh(m) {}
94  };
95 
96 
97  /* ********************************************************************* */
98  /* */
99  /* II. Interpolation of functions. */
100  /* */
101  /* ********************************************************************* */
102 
103 
104  template <typename VECT, typename F, typename M>
105  inline void interpolation_function__(const mesh_fem &mf, VECT &V,
106  F &f, const dal::bit_vector &dofs,
107  const M &, gmm::abstract_null_type) {
108  size_type Q = mf.get_qdim();
109  GMM_ASSERT1(gmm::vect_size(V) == mf.nb_basic_dof() && Q == 1,
110  "Dof vector has not the right size");
111  for (dal::bv_visitor i(dofs); !i.finished(); ++i)
112  V[i] = f(mf.point_of_basic_dof(i));
113  }
114 
115  template <typename VECT, typename F, typename M>
116  inline void interpolation_function__(const mesh_fem &mf, VECT &V,
117  F &f, const dal::bit_vector &dofs,
118  const M &v, gmm::abstract_vector) {
119  size_type N = gmm::vect_size(v), Q = mf.get_qdim();
120  GMM_ASSERT1(gmm::vect_size(V) == mf.nb_basic_dof()*N/Q,
121  "Dof vector has not the right size");
122  for (dal::bv_visitor i(dofs); !i.finished(); ++i)
123  if (i % Q == 0)
124  gmm::copy(f(mf.point_of_basic_dof(i)),
125  gmm::sub_vector(V, gmm::sub_interval(i*N/Q, N)));
126  }
127 
128  template <typename VECT, typename F, typename M>
129  inline void interpolation_function__(const mesh_fem &mf, VECT &V,
130  F &f, const dal::bit_vector &dofs,
131  const M &mm, gmm::abstract_matrix) {
132  // typedef typename gmm::linalg_traits<VECT>::value_type T;
133  size_type Nr = gmm::mat_nrows(mm), Nc = gmm::mat_ncols(mm), N = Nr*Nc;
134  size_type Q = mf.get_qdim();
135  base_matrix m(Nr, Nc);
136  GMM_ASSERT1(gmm::vect_size(V) == mf.nb_basic_dof()*N/Q,
137  "Dof vector has not the right size");
138  for (dal::bv_visitor i(dofs); !i.finished(); ++i)
139  if (i % Q == 0) {
140  gmm::copy(f(mf.point_of_basic_dof(i)), m);
141  for (size_type j = 0; j < Nc; ++j)
142  gmm::copy(gmm::mat_col(m, j),
143  gmm::sub_vector(V, gmm::sub_interval(i*N/Q+j*Nr, Nr)));
144  }
145 
146  }
147 
148  template <typename VECT, typename F, typename M>
149  inline void interpolation_function_(const mesh_fem &mf, VECT &V,
150  F &f, const dal::bit_vector &dofs,
151  const M &m) {
152  interpolation_function__(mf, V, f, dofs, m,
153  typename gmm::linalg_traits<M>::linalg_type());
154  }
155 
156 #if GETFEM_PARA_LEVEL > 0
157  template <typename T>
158  void take_one_op(void *a, void *b, int *len, MPI_Datatype *) {
159  T aa = *((T*)a);
160  return aa ? aa : *((T*)b);
161  }
162 
163  template <typename T>
164  inline MPI_Op mpi_take_one_op(T) {
165  static bool isinit = false;
166  static MPI_Op op;
167  if (!isinit) {
168  MPI_Op_create(take_one_op<T>, true, &op);
169  isinit = true;
170  }
171  return op;
172  }
173 #endif
174 
175  // TODO : verify that rhs is a lagrange fem
176  /**
177  @brief interpolation of a function f on mf_target.
178  - mf_target must be of lagrange type.
179  - mf_target's qdim should be equal to the size of the return value of f,
180  or equal to 1
181  - V should have the right size
182  CAUTION: with the parallized version (GETFEM_PARA_LEVEL >= 2) the
183  resulting vector V is distributed.
184  */
185  template <typename VECT, typename F>
186  void interpolation_function(mesh_fem &mf_target, const VECT &VV, F &f,
188  typedef typename gmm::linalg_traits<VECT>::value_type T;
189  size_type qqdimt = gmm::vect_size(VV) / mf_target.nb_dof();
190  std::vector<T> V(mf_target.nb_basic_dof()*qqdimt);
191  mf_target.linked_mesh().intersect_with_mpi_region(rg);
192  dal::bit_vector dofs = mf_target.basic_dof_on_region(rg);
193  if (dofs.card() > 0)
194  interpolation_function_(mf_target, V, f, dofs,
195  f(mf_target.point_of_basic_dof(dofs.first())));
196 
197  if (mf_target.is_reduced()) {
198  for (size_type k = 0; k < qqdimt; ++k)
199  gmm::mult(mf_target.reduction_matrix(),
200  gmm::sub_vector(V,
201  gmm::sub_slice(k, mf_target.nb_basic_dof(),
202  qqdimt)),
203  gmm::sub_vector(const_cast<VECT &>(VV),
204  gmm::sub_slice(k, mf_target.nb_dof(),
205  qqdimt)));
206  }
207  else
208  gmm::copy(V, const_cast<VECT &>(VV));
209  }
210 
211  /* ********************************************************************* */
212  /* */
213  /* III. Interpolation between two meshes. */
214  /* */
215  /* ********************************************************************* */
216 
217  /* ------------------------------ Interface -----------------------------*/
218 
219  /**
220  @brief interpolation/extrapolation of (mf_source, U) on mf_target.
221  - mf_target must be of lagrange type.
222  - mf_target's qdim should be equal to mf_source qdim, or equal to 1
223  - U.size() >= mf_source.get_qdim()
224  - V.size() >= (mf_target.nb_dof() / mf_target.get_qdim())
225  * mf_source.get_qdim()
226 
227  With extrapolation = 0 a strict interpolation is done, with extrapolation = 1
228  an extrapolation of the exterior points near the boundary is done (if any)
229  and with extrapolation = 2 all exterior points are extrapolated (could be expensive).
230 
231  If both mesh_fem shared the same mesh object, a fast interpolation
232  will be used.
233 
234  If rg_source and rg_target are provided the operation is restricted to
235  these regions. rg_source must contain only convexes.
236  */
237  template<typename VECTU, typename VECTV>
238  void interpolation(const mesh_fem &mf_source, const mesh_fem &mf_target,
239  const VECTU &U, VECTV &V, int extrapolation = 0,
240  double EPS = 1E-10,
241  mesh_region rg_source=mesh_region::all_convexes(),
242  mesh_region rg_target=mesh_region::all_convexes());
243 
244  /**
245  @brief Build the interpolation matrix of mf_source on mf_target.
246  the matrix M is
247  such that (V = M*U) == interpolation(mf_source, mf_target, U, V).
248 
249  Useful for repeated interpolations.
250  For performance reasons the matrix M is recommended to be either
251  a row or a row and column matrix.
252 
253  If rg_source and rg_target are provided the operation is restricted to
254  these regions. rg_source must contain only convexes.
255  */
256  template<typename MAT>
257  void interpolation(const mesh_fem &mf_source, const mesh_fem &mf_target,
258  MAT &M, int extrapolation = 0, double EPS = 1E-10,
259  mesh_region rg_source=mesh_region::all_convexes(),
260  mesh_region rg_target=mesh_region::all_convexes());
261 
262 
263  /* --------------------------- Implementation ---------------------------*/
264 
265  /*
266  interpolation of a solution on same mesh.
267  - &mf_target.linked_mesh() == &mf_source.linked_mesh()
268  - mf_target must be of lagrange type.
269  - mf_target's qdim should be equal to mf_source qdim, or equal to 1
270  - U.size() >= mf_source.get_qdim()
271  - V.size() >= (mf_target.nb_dof() / mf_target.get_qdim())
272  * mf_source.get_qdim()
273  */
274  template<typename VECTU, typename VECTV, typename MAT>
275  void interpolation_same_mesh(const mesh_fem &mf_source,
276  const mesh_fem &mf_target,
277  const VECTU &UU, VECTV &VV,
278  MAT &MM, int version) {
279 
280  typedef typename gmm::linalg_traits<VECTU>::value_type T;
281  base_matrix G;
282  dim_type qdim = mf_source.get_qdim();
283  dim_type qqdim = dim_type(gmm::vect_size(UU)/mf_source.nb_dof());
284  std::vector<T> val(qdim);
285  std::vector<std::vector<T> > coeff;
286  std::vector<size_type> dof_source;
287  GMM_ASSERT1(qdim == mf_target.get_qdim() || mf_target.get_qdim() == 1,
288  "Attempt to interpolate a field of dimension "
289  << qdim << " on a mesh_fem whose Qdim is " <<
290  int(mf_target.get_qdim()));
291  size_type qmult = mf_source.get_qdim()/mf_target.get_qdim();
292  size_type qqdimt = qqdim * mf_source.get_qdim()/mf_target.get_qdim();
293  fem_precomp_pool fppool;
294  std::vector<size_type> dof_t_passes(mf_target.nb_basic_dof());
295  std::vector<T> U(mf_source.nb_basic_dof()*qqdim);
296  std::vector<T> V(mf_target.nb_basic_dof()*qqdimt);
297  gmm::row_matrix<gmm::rsvector<scalar_type> >
298  M(mf_target.nb_basic_dof(), mf_source.nb_basic_dof());
299 
300  if (version == 0) mf_source.extend_vector(UU, U);
301 
302  /* we should sort convexes by their fem */
303  for (dal::bv_visitor cv(mf_source.convex_index()); !cv.finished(); ++cv) {
304  bgeot::pgeometric_trans pgt=mf_source.linked_mesh().trans_of_convex(cv);
305  pfem pf_s = mf_source.fem_of_element(cv);
306  if (!mf_target.convex_index().is_in(cv))
307  continue;
308  pfem pf_t = mf_target.fem_of_element(cv);
309  size_type nbd_s = pf_s->nb_dof(cv);
310  size_type nbd_t = pf_t->nb_dof(cv);
311  mesh_fem::ind_dof_ct::const_iterator itdof;
312  size_type cvnbdof = mf_source.nb_basic_dof_of_element(cv);
313 
314  bool discontinuous_source = false;
315  for (size_type dof=0; dof < nbd_s; ++dof)
316  if (!dof_linkable(pf_s->dof_types()[dof])) {
317  discontinuous_source = true;
318  break;
319  }
320 
321  if (version == 0) {
322  coeff.resize(qqdim);
323  for (size_type qq=0; qq < qqdim; ++qq) {
324  coeff[qq].resize(cvnbdof);
325  itdof = mf_source.ind_basic_dof_of_element(cv).begin();
326  for (size_type k = 0; k < cvnbdof; ++k, ++itdof) {
327  coeff[qq][k] = U[(*itdof)*qqdim+qq];
328  }
329  }
330  }
331  if (pf_s->need_G())
332  bgeot::vectors_to_base_matrix
333  (G, mf_source.linked_mesh().points_of_convex(cv));
334 
335  GMM_ASSERT1(pf_t->target_dim() == 1,
336  "won't interpolate on a vector FEM... ");
337  pfem_precomp pfp = fppool(pf_s, pf_t->node_tab(cv));
338  fem_interpolation_context ctx(pgt,pfp,size_type(-1), G, cv,
339  short_type(-1));
340  itdof = mf_target.ind_basic_dof_of_element(cv).begin();
341  if (version != 0) {
342  const mesh_fem::ind_dof_ct &idct
343  = mf_source.ind_basic_dof_of_element(cv);
344  dof_source.assign(idct.begin(), idct.end());
345  }
346  for (size_type i = 0; i < nbd_t; ++i, itdof+=mf_target.get_qdim()) {
347  size_type dof_t = *itdof*qmult;
348  if (!discontinuous_source && dof_t_passes[*itdof] > 0) continue;
349  dof_t_passes[*itdof] += 1;
350  ctx.set_ii(i);
351  if (version == 0) {
352  for (size_type qq=0; qq < qqdim; ++qq) {
353  pf_s->interpolation(ctx, coeff[qq], val, qdim);
354  for (size_type k=0; k < qdim; ++k)
355  V[(dof_t + k)*qqdim+qq] += val[k];
356  }
357  }
358  else {
359  base_matrix Mloc(qdim, mf_source.nb_basic_dof_of_element(cv));
360  pf_s->interpolation(ctx, Mloc, qdim);
361  for (size_type k=0; k < qdim; ++k) {
362  for (size_type j=0; j < dof_source.size(); ++j) {
363  M(dof_t + k, dof_source[j]) += Mloc(k, j);
364  }
365  }
366  }
367  }
368  }
369 
370  // calculate averages for discontinuous source and continuous target
371  for (size_type i = 0; i < mf_target.nb_basic_dof(); ++i) {
372  size_type dof_t = i*qmult;
373  scalar_type passes = scalar_type(dof_t_passes[i]);
374  if (version == 0 && passes > scalar_type(0))
375  for (size_type qq=0; qq < qqdim; ++qq)
376  for (size_type k=0; k < qdim; ++k)
377  V[(dof_t + k)*qqdim+qq] /= passes;
378  else if (passes > scalar_type(0))
379  for (size_type k=0; k < qdim; ++k)
380  for (size_type j=0; j < dof_source.size(); ++j)
381  gmm::scale(gmm::mat_row(M, dof_t + k), scalar_type(1)/passes);
382  }
383 
384  if (version == 0)
385  mf_target.reduce_vector(V, VV);
386  else {
387  if (mf_target.is_reduced())
388  if (mf_source.is_reduced()) {
389  gmm::row_matrix<gmm::rsvector<scalar_type> >
390  MMM(mf_target.nb_dof(), mf_source.nb_basic_dof());
391  gmm::mult(mf_target.reduction_matrix(), M, MMM);
392  gmm::mult(MMM, mf_source.extension_matrix(), MM);
393  }
394  else
395  gmm::mult(mf_target.reduction_matrix(), M, MM);
396  else
397  if (mf_source.is_reduced())
398  gmm::mult(M, mf_source.extension_matrix(), MM);
399  else
400  gmm::copy(M, MM);
401  }
402  }
403 
404 
405  /*
406  interpolation of a solution on another mesh.
407  - mti contains the points where to interpolate.
408  - the solution should be continuous.
409  */
410  template<typename VECTU, typename VECTV, typename MAT>
411  void interpolation(const mesh_fem &mf_source,
412  mesh_trans_inv &mti,
413  const VECTU &UU, VECTV &V, MAT &MM,
414  int version, int extrapolation = 0,
415  dal::bit_vector *dof_untouched = 0,
416  mesh_region rg_source=mesh_region::all_convexes()) {
417 
418  typedef typename gmm::linalg_traits<VECTU>::value_type T;
419  const mesh &msh(mf_source.linked_mesh());
420  dim_type qdim_s = mf_source.get_qdim();
421  dim_type qqdim = dim_type(gmm::vect_size(UU)/mf_source.nb_dof());
422 
423  std::vector<T> U(mf_source.nb_basic_dof()*qqdim);
424  gmm::row_matrix<gmm::rsvector<scalar_type> > M;
425  if (version != 0) M.resize(gmm::mat_nrows(MM), mf_source.nb_basic_dof());
426 
427  if (version == 0) mf_source.extend_vector(UU, U);
428 
429  mti.distribute(extrapolation, rg_source);
430  std::vector<size_type> itab;
431  base_matrix G;
432 
433  /* interpolation */
434  dal::bit_vector points_to_do; points_to_do.add(0, mti.nb_points());
435  std::vector<T> val(qdim_s);
436  std::vector<std::vector<T> > coeff;
437  base_tensor Z;
438  std::vector<size_type> dof_source;
439 
440  for (dal::bv_visitor cv(mf_source.convex_index()); !cv.finished(); ++cv) {
441  bgeot::pgeometric_trans pgt = msh.trans_of_convex(cv);
442  mti.points_on_convex(cv, itab);
443  if (itab.size() == 0) continue;
444 
445  pfem pf_s = mf_source.fem_of_element(cv);
446  if (pf_s->need_G())
447  bgeot::vectors_to_base_matrix(G, msh.points_of_convex(cv));
448 
449  fem_interpolation_context ctx(pgt, pf_s, base_node(), G, cv,
450  short_type(-1));
451  if (version == 0) {
452  coeff.resize(qqdim);
453  size_type cvnbdof = mf_source.nb_basic_dof_of_element(cv);
454  mesh_fem::ind_dof_ct::const_iterator itdof;
455  for (size_type qq=0; qq < qqdim; ++qq) {
456  coeff[qq].resize(cvnbdof);
457  itdof = mf_source.ind_basic_dof_of_element(cv).begin();
458  for (size_type k = 0; k < cvnbdof; ++k, ++itdof) {
459  coeff[qq][k] = U[(*itdof)*qqdim+qq];
460  }
461  }
462  }
463  if (version != 0) {
464  const mesh_fem::ind_dof_ct &idct
465  = mf_source.ind_basic_dof_of_element(cv);
466  dof_source.assign(idct.begin(), idct.end());
467  }
468  for (size_type i = 0; i < itab.size(); ++i) {
469  size_type ipt = itab[i];
470  if (points_to_do.is_in(ipt)) {
471  points_to_do.sup(ipt);
472  ctx.set_xref(mti.reference_coords()[ipt]);
473  size_type dof_t = mti.id_of_point(ipt);
474  size_type pos = dof_t * qdim_s;
475  if (version == 0) {
476  for (size_type qq=0; qq < qqdim; ++qq) {
477  pf_s->interpolation(ctx, coeff[qq], val, qdim_s);
478  for (size_type k=0; k < qdim_s; ++k)
479  V[(pos + k)*qqdim+qq] = val[k];
480  }
481  // Part to be improved if one wants in option to be able to
482  // interpolate the gradient.
483  // if (PVGRAD) {
484  // base_matrix grad(mdim, qdim);
485  // pf_s->interpolation_grad(ctx,coeff,gmm::transposed(grad), qdim);
486  // std::copy(grad.begin(), grad.end(), V.begin()+dof_t*qdim*mdim);
487  // }
488  } else {
489  base_matrix Mloc(qdim_s, mf_source.nb_basic_dof_of_element(cv));
490  pf_s->interpolation(ctx, Mloc, qdim_s);
491  for (size_type k=0; k < qdim_s; ++k) {
492  for (size_type j=0; j < gmm::mat_ncols(Mloc); ++j)
493  M(pos+k, dof_source[j]) = Mloc(k,j);
494  /* does not work with col matrices
495  gmm::clear(gmm::mat_row(M, pos+k));
496  gmm::copy(gmm::mat_row(Mloc, k),
497  gmm::sub_vector(gmm::mat_row(M, pos+k), isrc));
498  */
499  }
500  }
501  }
502  }
503  }
504  if (points_to_do.card() != 0) {
505  if (dof_untouched) {
506  dof_untouched->clear();
507  for (dal::bv_visitor ipt(points_to_do); !ipt.finished(); ++ipt)
508  dof_untouched->add(mti.id_of_point(ipt));
509  }
510  else {
511  dal::bit_vector dofs_to_do;
512  for (dal::bv_visitor ipt(points_to_do); !ipt.finished(); ++ipt)
513  dofs_to_do.add(mti.id_of_point(ipt));
514  GMM_WARNING2("in interpolation (different meshes),"
515  << dofs_to_do.card() << " dof of target mesh_fem have "
516  << " been missed\nmissing dofs : " << dofs_to_do);
517  }
518  }
519 
520  if (version != 0) {
521  if (mf_source.is_reduced())
522  gmm::mult(M, mf_source.extension_matrix(), MM);
523  else
524  gmm::copy(M, MM);
525  }
526 
527  }
528 
529  template<typename VECTU, typename VECTV>
530  void interpolation(const mesh_fem &mf_source, mesh_trans_inv &mti,
531  const VECTU &U, VECTV &V, int extrapolation = 0,
532  dal::bit_vector *dof_untouched = 0,
533  mesh_region rg_source=mesh_region::all_convexes()) {
534  base_matrix M;
535  GMM_ASSERT1((gmm::vect_size(U) % mf_source.nb_dof()) == 0 &&
536  gmm::vect_size(V)!=0, "Dimension of vector mismatch");
537  interpolation(mf_source, mti, U, V, M, 0, extrapolation, dof_untouched, rg_source);
538  }
539 
540 
541 
542  /*
543  interpolation of a solution on another mesh.
544  - mf_target must be of lagrange type.
545  - the solution should be continuous..
546  */
547  template<typename VECTU, typename VECTV, typename MAT>
548  void interpolation(const mesh_fem &mf_source, const mesh_fem &mf_target,
549  const VECTU &U, VECTV &VV, MAT &MM,
550  int version, int extrapolation,
551  double EPS,
552  mesh_region rg_source=mesh_region::all_convexes(),
553  mesh_region rg_target=mesh_region::all_convexes()) {
554 
555  //Check if it is a torus mesh_fem
556  const torus_mesh_fem * pmf_torus = dynamic_cast<const torus_mesh_fem *>(&mf_target);
557  if(pmf_torus != 0){
558  interpolation_to_torus_mesh_fem(mf_source, *pmf_torus,
559  U, VV, MM, version, extrapolation, EPS, rg_source, rg_target);
560  return;
561  }
562 
563  typedef typename gmm::linalg_traits<VECTU>::value_type T;
564  dim_type qqdim = dim_type(gmm::vect_size(U)/mf_source.nb_dof());
565  size_type qqdimt = qqdim * mf_source.get_qdim()/mf_target.get_qdim();
566  std::vector<T> V(mf_target.nb_basic_dof()*qqdimt);
567  mf_target.extend_vector(VV,V);
568  gmm::row_matrix<gmm::rsvector<scalar_type> > M;
569  if (version != 0) M.resize(mf_target.nb_basic_dof(), mf_source.nb_dof());
570 
571  const mesh &msh(mf_source.linked_mesh());
572  getfem::mesh_trans_inv mti(msh, EPS);
573  size_type qdim_s = mf_source.get_qdim(), qdim_t = mf_target.get_qdim();
574  GMM_ASSERT1(qdim_s == qdim_t || qdim_t == 1,
575  "Attempt to interpolate a field of dimension "
576  << qdim_s << " on a mesh_fem whose Qdim is " << qdim_t);
577 
578  /* test if the target mesh_fem is really of Lagrange type. */
579  for (dal::bv_visitor cv(mf_target.convex_index()); !cv.finished();++cv) {
580  pfem pf_t = mf_target.fem_of_element(cv);
581  GMM_ASSERT1(pf_t->target_dim() == 1 && pf_t->is_lagrange(),
582  "Target fem not convenient for interpolation");
583  }
584  /* initialisation of the mesh_trans_inv */
585  bool is_target_torus = dynamic_cast<const torus_mesh *>(&mf_target.linked_mesh());
586  if (rg_target.id() == mesh_region::all_convexes().id()) {
587  size_type nbpts = mf_target.nb_basic_dof() / qdim_t;
588  for (size_type i = 0; i < nbpts; ++i){
589  if (is_target_torus){
590  auto p = mf_target.point_of_basic_dof(i * qdim_t);
591  p.resize(msh.dim());
592  mti.add_point(p);
593  }
594  else mti.add_point(mf_target.point_of_basic_dof(i * qdim_t));
595  }
596  }
597  else {
598  for (dal::bv_visitor_c dof(mf_target.basic_dof_on_region(rg_target));
599  !dof.finished(); ++dof)
600  if (dof % qdim_t == 0){
601  if (is_target_torus){
602  auto p = mf_target.point_of_basic_dof(dof);
603  p.resize(msh.dim());
604  mti.add_point_with_id(p, dof/qdim_t);
605  }
606  else
607  mti.add_point_with_id(mf_target.point_of_basic_dof(dof),dof/qdim_t);
608  }
609  }
610  interpolation(mf_source, mti, U, V, M, version, extrapolation, 0,rg_source);
611 
612  if (version == 0)
613  mf_target.reduce_vector(V, VV);
614  else {
615  if (mf_target.is_reduced())
616  gmm::mult(mf_target.reduction_matrix(), M, MM);
617  else
618  gmm::copy(M, MM);
619  }
620 
621  }
622 
623  /*
624  interpolation of a solution on another torus mesh.
625  - the solution should be continuous.
626  */
627  template<typename VECTU, typename VECTV, typename MAT>
628  void interpolation_to_torus_mesh_fem(const mesh_fem &mf_source, const torus_mesh_fem &mf_target,
629  const VECTU &U, VECTV &VV, MAT &MM,
630  int version, int extrapolation,
631  double EPS,
632  mesh_region rg_source=mesh_region::all_convexes(),
633  mesh_region rg_target=mesh_region::all_convexes()) {
634 
635  typedef typename gmm::linalg_traits<VECTU>::value_type T;
636  dim_type qqdim = dim_type(gmm::vect_size(U)/mf_source.nb_dof());
637  size_type qqdimt = qqdim * mf_source.get_qdim()/mf_target.get_qdim();
638  std::vector<T> V(mf_target.nb_basic_dof()*qqdimt);
639  mf_target.extend_vector(VV,V);
640  gmm::row_matrix<gmm::rsvector<scalar_type> >
641  M(mf_target.nb_basic_dof(), mf_source.nb_dof());
642 
643  const mesh &msh(mf_source.linked_mesh());
644  getfem::mesh_trans_inv mti(msh, EPS);
645  size_type qdim_s = mf_source.get_qdim(), qdim_t = mf_target.get_qdim();
646  GMM_ASSERT1(qdim_s == qdim_t || qdim_t == 1,
647  "Attempt to interpolate a field of dimension "
648  << qdim_s << " on a mesh_fem whose Qdim is " << qdim_t);
649 
650  /* test if the target mesh_fem is convenient for interpolation. */
651  for (dal::bv_visitor cv(mf_target.convex_index()); !cv.finished();++cv) {
652  pfem pf_t = mf_target.fem_of_element(cv);
653 
654  GMM_ASSERT1(pf_t->target_dim() == 1 ||
655  (mf_target.get_qdim() == mf_target.linked_mesh().dim()),
656  "Target fem not convenient for interpolation");
657  }
658  /* initialisation of the mesh_trans_inv */
659  if (rg_target.id() == mesh_region::all_convexes().id()) {
660  size_type nbpts = mf_target.nb_basic_dof() / qdim_t;
661  for (size_type i = 0; i < nbpts; ++i)
662  {
663  base_node p(msh.dim());
664  for (size_type k=0; k < msh.dim(); ++k) p[k] = mf_target.point_of_basic_dof(i * qdim_t)[k];
665  mti.add_point(p);
666  }
667  interpolation(mf_source, mti, U, V, M, version, extrapolation);
668  }
669  else {
670  for (dal::bv_visitor_c dof(mf_target.basic_dof_on_region(rg_target)); !dof.finished(); ++dof)
671  if (dof % qdim_t == 0)
672  {
673  base_node p(msh.dim());
674  for (size_type k=0; k < msh.dim(); ++k) p[k] = mf_target.point_of_basic_dof(dof)[k];
675  mti.add_point_with_id(p, dof/qdim_t);
676  }
677  interpolation(mf_source, mti, U, V, M, version, extrapolation, 0, rg_source);
678  }
679 
680  if (version == 0)
681  mf_target.reduce_vector(V, VV);
682  else {
683  if (mf_target.is_reduced())
684  gmm::mult(mf_target.reduction_matrix(), M, MM);
685  else
686  gmm::copy(M, MM);
687  }
688  }
689 
690 
691 
692  template<typename VECTU, typename VECTV>
693  void interpolation(const mesh_fem &mf_source, const mesh_fem &mf_target,
694  const VECTU &U, VECTV &V, int extrapolation,
695  double EPS,
696  mesh_region rg_source, mesh_region rg_target) {
697  base_matrix M;
698  GMM_ASSERT1((gmm::vect_size(U) % mf_source.nb_dof()) == 0
699  && (gmm::vect_size(V) % mf_target.nb_dof()) == 0
700  && gmm::vect_size(V) != 0, "Dimensions mismatch");
701  if (&mf_source.linked_mesh() == &mf_target.linked_mesh() &&
702  rg_source.id() == mesh_region::all_convexes().id() &&
703  rg_target.id() == mesh_region::all_convexes().id())
704  interpolation_same_mesh(mf_source, mf_target, U, V, M, 0);
705  else {
706  omp_distribute<VECTV> V_distributed;
707  auto partitioning_allowed = rg_source.is_partitioning_allowed();
708  rg_source.prohibit_partitioning();
710  auto &V_thrd = V_distributed.thrd_cast();
711  gmm::resize(V_thrd, V.size());
713  mf_source, mf_target, U, V_thrd, M, 0, extrapolation, EPS,
714  rg_source, rg_target);
715  )
716  for (size_type thread=0; thread != V_distributed.num_threads(); ++thread){
717  auto &V_thrd2 = V_distributed(thread);
718  for (size_type i = 0; i < V_thrd2.size(); ++i) {
719  if (gmm::abs(V_thrd2[i]) > EPS) V[i] = V_thrd2[i];
720  }
721  }
722  if (partitioning_allowed) rg_source.allow_partitioning();
723  }
724  }
725 
726  template<typename MAT>
727  void interpolation(const mesh_fem &mf_source, const mesh_fem &mf_target,
728  MAT &M, int extrapolation, double EPS,
729  mesh_region rg_source, mesh_region rg_target) {
730  GMM_ASSERT1(mf_source.nb_dof() == gmm::mat_ncols(M)
731  && (gmm::mat_nrows(M) % mf_target.nb_dof()) == 0
732  && gmm::mat_nrows(M) != 0, "Dimensions mismatch");
733  std::vector<scalar_type> U, V;
734  if (&mf_source.linked_mesh() == &mf_target.linked_mesh() &&
735  rg_source.id() == mesh_region::all_convexes().id() &&
736  rg_target.id() ==mesh_region::all_convexes().id())
737  interpolation_same_mesh(mf_source, mf_target, U, V, M, 1);
738  else
739  interpolation(mf_source, mf_target, U, V, M, 1, extrapolation, EPS,
740  rg_source, rg_target);
741  }
742 
743 
744  /**Interpolate mesh_fem data to im_data.
745  The qdim of mesh_fem must be equal to im_data nb_tensor_elem.
746  Both im_data and mesh_fem must reside in the same mesh.
747  Only convexes defined with both mesh_fem and im_data will be interpolated.
748  The use_im_data_filter flag controls the use of the filtered region of
749  im_data (default) or the use the full mesh.
750  */
751  template <typename VECT>
752  void interpolation_to_im_data(const mesh_fem &mf_source, const im_data &im_target,
753  const VECT &nodal_data, VECT &int_pt_data,
754  bool use_im_data_filter = true) {
755  // typedef typename gmm::linalg_traits<const VECT>::value_type T;
756 
757  dim_type qdim = mf_source.get_qdim();
758  size_type nb_dof = mf_source.nb_dof();
759  size_type nb_basic_dof = mf_source.nb_basic_dof();
760  size_type nodal_data_size = gmm::vect_size(nodal_data);
761  dim_type data_qdim = nodal_data_size / nb_dof;
762 
763  GMM_ASSERT1(data_qdim * mf_source.nb_dof() == nodal_data_size,
764  "Incompatible size of mesh fem " << mf_source.nb_dof() * data_qdim <<
765  " with the data " << nodal_data_size);
766 
767  GMM_ASSERT1(qdim * data_qdim == im_target.nb_tensor_elem(),
768  "Incompatible size of qdim for mesh_fem " << qdim
769  << " and im_data " << im_target.nb_tensor_elem());
770  GMM_ASSERT1(&mf_source.linked_mesh() == &im_target.linked_mesh(),
771  "mf_source and im_data do not share the same mesh.");
772 
773  GMM_ASSERT1(nb_dof * data_qdim == nodal_data_size,
774  "Provided nodal data size is " << nodal_data_size
775  << " but expecting vector size of " << nb_dof);
776 
777  size_type size_im_data = im_target.nb_index(use_im_data_filter)
778  * im_target.nb_tensor_elem();
779  GMM_ASSERT1(size_im_data == gmm::vect_size(int_pt_data),
780  "Provided im data size is " << gmm::vect_size(int_pt_data)
781  << " but expecting vector size of " << size_im_data);
782 
783  VECT extended_nodal_data_((nb_dof != nb_basic_dof) ? nb_basic_dof * data_qdim : 0);
784  if (nb_dof != nb_basic_dof)
785  mf_source.extend_vector(nodal_data, extended_nodal_data_);
786  const VECT &extended_nodal_data = (nb_dof == nb_basic_dof) ? nodal_data : extended_nodal_data_;
787 
788  dal::bit_vector im_data_convex_index(im_target.convex_index(use_im_data_filter));
789 
790  base_matrix G;
791  base_vector coeff;
792  bgeot::base_tensor tensor_int_point(im_target.tensor_size());
793  fem_precomp_pool fppool;
794  for (dal::bv_visitor cv(im_data_convex_index); !cv.finished(); ++cv) {
795 
796  bgeot::pgeometric_trans pgt = mf_source.linked_mesh().trans_of_convex(cv);
797  pfem pf_source = mf_source.fem_of_element(cv);
798  if (pf_source == NULL)
799  continue;
800 
801  mesh_fem::ind_dof_ct::const_iterator it_dof;
802  size_type cv_nb_dof = mf_source.nb_basic_dof_of_element(cv);
803  size_type nb_nodal_pt = cv_nb_dof / qdim;
804  coeff.resize(cv_nb_dof);
805  getfem::slice_vector_on_basic_dof_of_element(mf_source, extended_nodal_data, cv, coeff);
806 
807  const getfem::papprox_integration pim(im_target.approx_int_method_of_element(cv));
808  if (pf_source->need_G())
809  bgeot::vectors_to_base_matrix(G, *(pim->pintegration_points()));
810 
811  pfem_precomp pfp = fppool(pf_source, pim->pintegration_points());
812 
813  // interior of the convex
814  size_type nb_int_pts;
815  size_type int_pt_id = im_target.index_of_first_point(cv, short_type(-1),
816  use_im_data_filter);
817  if (int_pt_id != size_type(-1)) {
818  nb_int_pts = im_target.nb_points_of_element(cv, short_type(-1));
819  fem_interpolation_context ctx(pgt, pfp, size_type(-1), G, cv);
820  for (size_type i = 0; i < nb_int_pts; ++i, ++int_pt_id) {
821  ctx.set_ii(i);
822  ctx.pf()->interpolation(ctx, coeff, tensor_int_point.as_vector(), qdim * data_qdim);
823  im_target.set_tensor(int_pt_data, int_pt_id, tensor_int_point);
824  }
825  }
826 
827  // convex faces
828  for (short_type f=0, nb_faces=im_target.nb_faces_of_element(cv);
829  f < nb_faces; ++f) {
830  int_pt_id = im_target.index_of_first_point(cv, f, use_im_data_filter);
831  if (int_pt_id != size_type(-1)) {
832  nb_int_pts = im_target.nb_points_of_element(cv, f);
833  fem_interpolation_context ctx(pgt, pfp, size_type(-1), G, cv, f);
834  size_type i0 = pim->ind_first_point_on_face(f);
835  for (size_type i = 0; i < nb_int_pts; ++i, ++int_pt_id) {
836  ctx.set_ii(i+i0);
837  ctx.pf()->interpolation(ctx, coeff, tensor_int_point.as_vector(), qdim * data_qdim);
838  im_target.set_tensor(int_pt_data, int_pt_id, tensor_int_point);
839  }
840  }
841  }
842  }//end of convex loop
843  }
844 
845 } /* end of namespace getfem. */
846 
847 
848 #endif /* GETFEM_INTERPOLATION_H__ */
getfem::im_data::linked_mesh
const mesh & linked_mesh() const
linked mesh
Definition: getfem_im_data.h:159
getfem::mesh_fem::get_qdim
virtual dim_type get_qdim() const
Return the Q dimension.
Definition: getfem_mesh_fem.h:312
bgeot::geotrans_inv
handles the geometric inversion for a given (supposedly quite large) set of points
Definition: bgeot_geotrans_inv.h:178
bgeot_torus.h
Provides mesh of torus.
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::im_data::nb_points_of_element
size_type nb_points_of_element(size_type cv, bool use_filter=false) const
Total number of points in element cv.
Definition: getfem_im_data.cc:65
getfem::fem_precomp_pool
handle a pool (i.e.
Definition: getfem_fem.h:717
getfem::im_data::nb_tensor_elem
size_type nb_tensor_elem() const
sum of tensor elements, M(3,3) will have 3*3=9 elements
Definition: getfem_im_data.cc:229
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
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_region::allow_partitioning
void allow_partitioning()
In multithreaded part of the program makes only a partition of the region visible in the index() and ...
Definition: getfem_mesh_region.cc:204
getfem::mesh_region::prohibit_partitioning
void prohibit_partitioning()
Disregard partitioning, which means being able to see the whole region in multithreaded code.
Definition: getfem_mesh_region.cc:220
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::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
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
dal_tree_sorted.h
a balanced tree stored in a dal::dynamic_array
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
getfem::fem_interpolation_context
structure passed as the argument of fem interpolation functions.
Definition: getfem_fem.h:749
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::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
getfem_mesh_fem.h
Define the getfem::mesh_fem class.
getfem::im_data::convex_index
dal::bit_vector convex_index(bool use_filter=false) const
List of convexes.
Definition: getfem_im_data.cc:159
getfem::omp_distribute
Use this template class for any object you want to distribute to open_MP threads.
Definition: getfem_omp.h:325
bgeot::geotrans_interpolation_context::set_ii
void set_ii(size_type ii__)
change the current point (assuming a geotrans_precomp_ is used)
Definition: bgeot_geometric_trans.h:463
getfem::im_data::index_of_first_point
size_type index_of_first_point(size_type cv, short_type f=short_type(-1), bool use_filter=false) const
Returns the index of the first integration point with no filtering.
Definition: getfem_im_data.cc:140
getfem_im_data.h
Provides indexing of integration points for mesh_im.
GETFEM_OMP_PARALLEL
#define GETFEM_OMP_PARALLEL(body)
Organizes a proper parallel omp section:
Definition: getfem_omp.h:483
getfem::dof_linkable
bool dof_linkable(pdof_description)
Says if the dof is linkable.
Definition: getfem_fem.cc:615
getfem::im_data::set_tensor
void set_tensor(VECT &V1, size_type cv, size_type i, const TENSOR &T, bool use_filter=true) const
set a tensor of an integration point from a raw vector data, described by the tensor size.
Definition: getfem_im_data.h:352
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
bgeot::geotrans_inv::add_point
size_type add_point(base_node p)
Add point p to the list of points.
Definition: bgeot_geotrans_inv.h:196
bgeot
Basic Geometric Tools.
Definition: bgeot_convex_ref.cc:27
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::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::im_data
im_data provides indexing to the integration points of a mesh im object.
Definition: getfem_im_data.h:69
bgeot::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:186
getfem::mesh_region
structure used to hold a set of convexes and/or convex faces.
Definition: getfem_mesh_region.h:55
getfem::im_data::nb_faces_of_element
short_type nb_faces_of_element(size_type cv) const
Number of (active) faces in element cv.
Definition: getfem_im_data.cc:104
getfem::fem_interpolation_context::pf
const pfem pf() const
get the current FEM descriptor
Definition: getfem_fem.h:782
getfem::interpolation_to_im_data
void interpolation_to_im_data(const mesh_fem &mf_source, const im_data &im_target, const VECT &nodal_data, VECT &int_pt_data, bool use_im_data_filter=true)
Interpolate mesh_fem data to im_data.
Definition: getfem_interpolation.h:752
getfem::interpolation_function
void interpolation_function(mesh_fem &mf_target, const VECT &VV, F &f, mesh_region rg=mesh_region::all_convexes())
interpolation of a function f on mf_target.
Definition: getfem_interpolation.h:186
getfem_torus.h
Provides mesh and mesh fem of torus.
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::im_data::nb_index
size_type nb_index(bool use_filter=false) const
Total numbers of index (integration points)
Definition: getfem_im_data.cc:57