Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
interpolate.hpp
Go to the documentation of this file.
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2008-01-31
7 
8  Copyright (C) 2008-2011 Universite Joseph Fourier (Grenoble I)
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 3.0 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public
21  License along with this library; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
29 #ifndef __interpolate_H
30 #define __interpolate_H 1
31 
32 namespace Feel
33 {
34 enum { INTERPOLATE_DIFFERENT_MESH=0, INTERPOLATE_SAME_MESH = 1 };
35 
49 template<typename SpaceType, typename FunctionType>
50 void
51 interpolate( boost::shared_ptr<SpaceType> const& space,
52  FunctionType const& f,
53  typename SpaceType::element_type& interp, int same_mesh = INTERPOLATE_DIFFERENT_MESH )
54 {
55  typedef typename SpaceType::value_type value_type;
56  typedef boost::multi_array<value_type,3> array_type;
57  typedef typename SpaceType::element_type interp_element_type;
58 
59  typedef typename SpaceType::mesh_type mesh_type;
60  typedef typename mesh_type::element_type geoelement_type;
61  typedef typename mesh_type::element_iterator mesh_element_iterator;
62 
63  typedef typename FunctionType::functionspace_type::mesh_type domain_mesh_type;
64  typedef typename domain_mesh_type::element_type domain_geoelement_type;
65  typedef typename domain_mesh_type::element_iterator domain_mesh_element_iterator;
66  // geometric mapping context
67  typedef typename mesh_type::gm_type gm_type;
68  typedef boost::shared_ptr<gm_type> gm_ptrtype;
69  typedef typename gm_type::template Context<vm::POINT|vm::GRAD|vm::KB|vm::JACOBIAN, geoelement_type> gmc_type;
70  typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
71 
72  typedef typename domain_mesh_type::gm_type domain_gm_type;
73  typedef boost::shared_ptr<domain_gm_type> domain_gm_ptrtype;
74  typedef typename domain_gm_type::template Context<vm::POINT|vm::GRAD|vm::KB|vm::JACOBIAN, domain_geoelement_type> domain_gmc_type;
75  typedef boost::shared_ptr<domain_gmc_type> domain_gmc_ptrtype;
76 
77  typedef typename FunctionType::functionspace_type::fe_type f_fe_type;
79  typedef boost::shared_ptr<f_fectx_type> f_fectx_ptrtype;
80 
81  // dof
82  typedef typename SpaceType::dof_type dof_type;
83 
84  // basis
85  typedef typename SpaceType::basis_type basis_type;
86 
87 
88  const bool same_basis = boost::is_same<basis_type, typename FunctionType::functionspace_type::basis_type>::value;
89  DVLOG(2) << "[interpolate] are the basis the same " << same_basis << "\n";
90  DVLOG(2) << "[interpolate] are the meshes the same " << same_mesh << "\n";
91 
92  // if same space type and mesh then return the function itself
93  if ( same_basis && same_mesh == INTERPOLATE_SAME_MESH )
94  {
95  DVLOG(2) << "[interpolate] Same mesh and same space\n";
96  interp = f;
97  return;
98  }
99 
100  dof_type const* __dof = space->dof().get();
101  basis_type const* __basis = space->basis().get();
102  gm_ptrtype __gm = space->gm();
103  typedef typename gm_type::precompute_ptrtype geopc_ptrtype;
104  typedef typename gm_type::precompute_type geopc_type;
105  geopc_ptrtype __geopc( new geopc_type( __gm, __basis->dual().points() ) );
106 
107 
108  f.updateGlobalValues();
109 
110  auto it = f.functionSpace()->mesh()->beginElementWithProcessId( space->mesh()->comm().rank() );
111  auto en = f.functionSpace()->mesh()->endElementWithProcessId( space->mesh()->comm().rank() );
112  if ( it==en ) return;
113 
114  //gmc_ptrtype __c( new gmc_type( __gm, *it, __geopc ) );
115 
116  //f.id( *fectx, fvalues );
117 
118  // if same mesh but not same function space (different order)
119  //if ( f.functionSpace()->mesh() == space->mesh() )
120  //if ( same_mesh == INTERPOLATE_SAME_MESH )
121  if ( ( MeshBase* )f.functionSpace()->mesh().get() == ( MeshBase* )space->mesh().get() )
122  {
123  DVLOG(2) << "[interpolate] Same mesh but not same space\n";
124 
125  domain_gm_ptrtype __dgm = f.functionSpace()->gm();
126  typedef typename domain_gm_type::precompute_ptrtype domain_geopc_ptrtype;
127  typedef typename domain_gm_type::precompute_type domain_geopc_type;
128  domain_geopc_ptrtype __dgeopc( new domain_geopc_type( __dgm, __basis->dual().points() ) );
129 
130  domain_gmc_ptrtype __c( new domain_gmc_type( __dgm, *it, __dgeopc ) );
131  auto pc = f.functionSpace()->fe()->preCompute( f.functionSpace()->fe(), __c->xRefs() );
132 
133  f_fectx_ptrtype fectx( new f_fectx_type( f.functionSpace()->fe(),
134  __c,
135  pc ) );
136 
137  typedef boost::multi_array<typename f_fectx_type::id_type,1> array_type;
138  array_type fvalues( f.idExtents( *fectx ) );
139 
140  for ( ; it != en; ++ it )
141  {
142  __c->update( *it );
143  fectx->update( __c, pc );
144  std::fill( fvalues.data(), fvalues.data()+fvalues.num_elements(), f_fectx_type::id_type::Zero() );
145  f.id( *fectx, fvalues );
146 
147  //std::cout << "interpfunc : " << interpfunc << "\n";
148  for ( uint16_type l = 0; l < basis_type::nLocalDof; ++l )
149  {
150 
151  const int ncdof = basis_type::is_product?basis_type::nComponents:1;
152 
153  for ( uint16_type comp = 0; comp < ncdof; ++comp )
154  {
155  size_type globaldof = boost::get<0>( __dof->localToGlobal( it->id(),
156  l, comp ) );
157 
158 #if 0
159  size_type globaldof_f = boost::get<0>( f.functionSpace()->dof()->localToGlobal( it->id(),l, 0 ) );
160  std::cout << "elt : " << it->id() << "\n"
161  << " l : " << l << "\n"
162  << " comp: " << comp << "\n"
163  << " dof: " << globaldof_f << "\n"
164  << " value: " << f( globaldof_f ) << "\n";
165 #endif
166 
167  //DVLOG(2) << "globaldof = " << globaldof << " firstldof = " << interp.firstLocalIndex() << " lastldof " << interp.lastLocalIndex() << "\n";
168  // update only values on the processor
169  if ( globaldof >= interp.firstLocalIndex() &&
170  globaldof < interp.lastLocalIndex() )
171  {
172  interp( globaldof ) = fvalues[l]( comp,0 );
173  //DVLOG(2) << "interp( " << globaldof << ")=" << interp( globaldof ) << "\n";
174  //std::cout << "interp( " << globaldof << ")=" << interp( globaldof ) << "\n";
175  }
176  }
177  }
178  }
179 
180  DVLOG(2) << "[interpolate] Same mesh but not same space done\n";
181  } // same mesh
182 
183  else // INTERPOLATE_DIFFERENT_MESH
184  {
185  DVLOG(2) << "[interpolate] different meshes\n";
186  domain_gm_ptrtype __dgm = f.functionSpace()->gm();
187  typedef typename domain_gm_type::precompute_ptrtype domain_geopc_ptrtype;
188  typedef typename domain_gm_type::precompute_type domain_geopc_type;
189  // get only one point
190  typename matrix_node<value_type>::type pts( mesh_type::nDim, 1 );
191  ublas::column( pts, 0 ) = ublas::column( __basis->dual().points(), 0 );
192 
193  domain_geopc_ptrtype __dgeopc( new domain_geopc_type( __dgm, pts ) );
194 
195  domain_gmc_ptrtype __c( new domain_gmc_type( __dgm, *it, __dgeopc ) );
196  auto pc = f.functionSpace()->fe()->preCompute( f.functionSpace()->fe(), __c->xRefs() );
197 
198  f_fectx_ptrtype fectx( new f_fectx_type( f.functionSpace()->fe(),
199  __c,
200  pc ) );
201  typedef boost::multi_array<typename f_fectx_type::id_type,1> array_type;
202  array_type fvalues( f.idExtents( *fectx ) );
203 
204 
205  typename domain_mesh_type::Inverse meshinv( f.functionSpace()->mesh() );
206 
207  /* initialisation of the mesh::inverse data structure */
208  typename SpaceType::dof_type::dof_points_const_iterator it_dofpt = space->dof()->dofPointBegin();
209  typename SpaceType::dof_type::dof_points_const_iterator en_dofpt = space->dof()->dofPointEnd();
210  size_type nbpts = 0;
211 
212  for ( ; it_dofpt != en_dofpt; ++it_dofpt, ++nbpts )
213  {
214  meshinv.addPointWithId( *it_dofpt );
215  }
216 
217  FEELPP_ASSERT( meshinv.nPoints() == nbpts )( meshinv.nPoints() )( nbpts ).error( "invalid number of points " );
218  meshinv.distribute();
219 
220  std::vector<bool> dof_done( nbpts );
221  std::fill( dof_done.begin(), dof_done.end(), false );
222  std::vector<boost::tuple<size_type,uint16_type > > itab;
223 
224  size_type first_dof = space->dof()->firstDof();
225 
226  for ( ; it != en; ++ it )
227  {
228  __c->update( *it );
229  meshinv.pointsInConvex( it->id(), itab );
230 
231  if ( itab.size() == 0 )
232  continue;
233 
234  for ( size_type i = 0; i < itab.size(); ++i )
235  {
236  // get dof id in target dof table
237  size_type dof;
238  uint16_type comp;
239  boost::tie( dof, comp ) = itab[i];
240 #if !defined( NDEBUG )
241  DVLOG(2) << "[interpolate] element : " << it->id() << " npts: " << itab.size() << " ptid: " << i
242  << " gdof: " << dof << " comp = " << comp << "\n";
243 #endif
244 
245  if ( !dof_done[dof-first_dof] )
246  {
247  dof_done[dof-first_dof]=true;
248  ublas::column( pts, 0 ) = meshinv.referenceCoords().find(dof)->second;
249  __dgeopc->update( pts );
250  //std::cout << "------------------------------------------------------------\n";
251  //std::cout << "pts = " << pts << "\n";
252  __c->update( *it );
253  pc->update( __c->xRefs() );
254  fectx->update( __c, pc );
255  //typename FunctionType::pc_type pc( f.functionSpace()->fe(), __c->xRefs() );
256  //typename FunctionType::id_type interpfunc( f.id( *__c, pc ) );
257  //typename FunctionType::id_type interpfunc;
258 
259  std::fill( fvalues.data(), fvalues.data()+fvalues.num_elements(), f_fectx_type::id_type::Zero() );
260  f.id( *fectx, fvalues );
261  //std::cout << "interpfunc : " << interpfunc << "\n";
262 
263  //for ( uint16_type comp = 0;comp < basis_type::nComponents;++comp )
264  {
265  //size_type globaldof = basis_type::nLocalDof*comp+first_dof+dof;
266  //size_type globaldof = first_dof+ndofcomp*comp+dof;
267  //size_type globaldof = first_dof+dof;
268  size_type globaldof = dof;
269 
270  // update only values on the processor
271  if ( globaldof >= interp.firstLocalIndex() &&
272  globaldof < interp.lastLocalIndex() )
273  interp( globaldof ) = fvalues[0]( comp,0 );
274 
275  //interp( globaldof ) = interpfunc(comp,0,0);
276  }
277  }
278  }
279  }
280 
281 #if 0
282 
283  for ( size_type i = 0; i < dof_done.size(); ++i )
284  {
285  if ( dof_done[i] != true )
286  {
287  LOG(INFO) << "[interpolate] dof not treated\n";
288  //FEELPP_ASSERT( dof_done[i] == true )( i ).warn ( "invalid dof, was not treated" );
289 
290  typename SpaceType::dof_type::dof_points_const_iterator it_dofpt = space->dof()->dofPointBegin();
291  typename SpaceType::dof_type::dof_points_const_iterator en_dofpt = space->dof()->dofPointEnd();
292  size_type nbpts = 0;
293 
294  for ( ; it_dofpt != en_dofpt; ++it_dofpt, ++nbpts )
295  {
296  //meshinv.addPointWithId( *it_dofpt );
297 
298  // be careful with indices in parallel
299  if ( boost::get<1>( *it_dofpt ) == i )
300  {
301  LOG(INFO) << " id : " << boost::get<1>( *it_dofpt ) << "\n";
302  LOG(INFO) << "coord : " << boost::get<0>( *it_dofpt ) << "\n";
303  LOG(INFO) << " comp : " << boost::get<2>( *it_dofpt ) << "\n";
304 
305  LOG(INFO) << "f( " << boost::get<0>( *it_dofpt ) << ")=" << f( boost::get<0>( *it_dofpt ) ) << "\n";
306  }
307  }
308  }
309  }
310 
311 #endif
312  }
313 
314  //std::cout << "interp=" << interp << "\n";
315 } // interpolate
316 
317 }
318 
319 #endif /* __interpolate_H */

Generated on Fri Oct 25 2013 14:24:16 for Feel++ by doxygen 1.8.4