29 #ifndef __interpolate_H
30 #define __interpolate_H 1
34 enum { INTERPOLATE_DIFFERENT_MESH=0, INTERPOLATE_SAME_MESH = 1 };
49 template<
typename SpaceType,
typename FunctionType>
52 FunctionType
const& f,
53 typename SpaceType::element_type& interp,
int same_mesh = INTERPOLATE_DIFFERENT_MESH )
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;
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;
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;
67 typedef typename mesh_type::gm_type gm_type;
68 typedef boost::shared_ptr<gm_type> gm_ptrtype;
70 typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
72 typedef typename domain_mesh_type::gm_type domain_gm_type;
73 typedef boost::shared_ptr<domain_gm_type> domain_gm_ptrtype;
75 typedef boost::shared_ptr<domain_gmc_type> domain_gmc_ptrtype;
77 typedef typename FunctionType::functionspace_type::fe_type f_fe_type;
79 typedef boost::shared_ptr<f_fectx_type> f_fectx_ptrtype;
82 typedef typename SpaceType::dof_type dof_type;
85 typedef typename SpaceType::basis_type basis_type;
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";
93 if ( same_basis && same_mesh == INTERPOLATE_SAME_MESH )
95 DVLOG(2) <<
"[interpolate] Same mesh and same space\n";
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() ) );
108 f.updateGlobalValues();
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;
121 if ( (
MeshBase* )f.functionSpace()->mesh().get() == (
MeshBase* )space->mesh().get() )
123 DVLOG(2) <<
"[interpolate] Same mesh but not same space\n";
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() ) );
130 domain_gmc_ptrtype __c(
new domain_gmc_type( __dgm, *it, __dgeopc ) );
131 auto pc = f.functionSpace()->fe()->preCompute( f.functionSpace()->fe(), __c->xRefs() );
133 f_fectx_ptrtype fectx(
new f_fectx_type( f.functionSpace()->fe(),
137 typedef boost::multi_array<typename f_fectx_type::id_type,1> array_type;
138 array_type fvalues( f.idExtents( *fectx ) );
140 for ( ; it != en; ++ 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 );
148 for ( uint16_type l = 0; l < basis_type::nLocalDof; ++l )
151 const int ncdof = basis_type::is_product?basis_type::nComponents:1;
153 for ( uint16_type comp = 0; comp < ncdof; ++comp )
155 size_type globaldof = boost::get<0>( __dof->localToGlobal( it->id(),
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";
169 if ( globaldof >= interp.firstLocalIndex() &&
170 globaldof < interp.lastLocalIndex() )
172 interp( globaldof ) = fvalues[l]( comp,0 );
180 DVLOG(2) <<
"[interpolate] Same mesh but not same space done\n";
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;
191 ublas::column( pts, 0 ) = ublas::column( __basis->dual().points(), 0 );
193 domain_geopc_ptrtype __dgeopc(
new domain_geopc_type( __dgm, pts ) );
195 domain_gmc_ptrtype __c(
new domain_gmc_type( __dgm, *it, __dgeopc ) );
196 auto pc = f.functionSpace()->fe()->preCompute( f.functionSpace()->fe(), __c->xRefs() );
198 f_fectx_ptrtype fectx(
new f_fectx_type( f.functionSpace()->fe(),
201 typedef boost::multi_array<typename f_fectx_type::id_type,1> array_type;
202 array_type fvalues( f.idExtents( *fectx ) );
205 typename domain_mesh_type::Inverse meshinv( f.functionSpace()->mesh() );
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();
212 for ( ; it_dofpt != en_dofpt; ++it_dofpt, ++nbpts )
214 meshinv.addPointWithId( *it_dofpt );
217 FEELPP_ASSERT( meshinv.nPoints() == nbpts )( meshinv.nPoints() )( nbpts ).error(
"invalid number of points " );
218 meshinv.distribute();
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;
224 size_type first_dof = space->dof()->firstDof();
226 for ( ; it != en; ++ it )
229 meshinv.pointsInConvex( it->id(), itab );
231 if ( itab.size() == 0 )
234 for (
size_type i = 0; i < itab.size(); ++i )
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";
245 if ( !dof_done[dof-first_dof] )
247 dof_done[dof-first_dof]=
true;
248 ublas::column( pts, 0 ) = meshinv.referenceCoords().find(dof)->second;
249 __dgeopc->update( pts );
253 pc->update( __c->xRefs() );
254 fectx->update( __c, pc );
259 std::fill( fvalues.data(), fvalues.data()+fvalues.num_elements(), f_fectx_type::id_type::Zero() );
260 f.id( *fectx, fvalues );
271 if ( globaldof >= interp.firstLocalIndex() &&
272 globaldof < interp.lastLocalIndex() )
273 interp( globaldof ) = fvalues[0]( comp,0 );
283 for (
size_type i = 0; i < dof_done.size(); ++i )
285 if ( dof_done[i] !=
true )
287 LOG(INFO) <<
"[interpolate] dof not treated\n";
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();
294 for ( ; it_dofpt != en_dofpt; ++it_dofpt, ++nbpts )
299 if ( boost::get<1>( *it_dofpt ) == i )
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";
305 LOG(INFO) <<
"f( " << boost::get<0>( *it_dofpt ) <<
")=" << f( boost::get<0>( *it_dofpt ) ) <<
"\n";