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
pointsetinterpolation.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: 2006-03-03
7 
8  Copyright (C) 2006 EPFL
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 __PointSetInterpolation_H
30 #define __PointSetInterpolation_H 1
31 
33 #include <boost/numeric/ublas/matrix.hpp>
34 #include <boost/numeric/ublas/io.hpp>
35 
36 namespace Feel
37 {
38 namespace ublas = boost::numeric::ublas;
39 
48 template<uint16_type Dim,
49  uint16_type Order,
50  typename T,
51  template<uint16_type,uint16_type,uint16_type> class Convex = Simplex>
52 class PointSetInterpolation : public PointSet<Convex<Dim,Order,Dim>,T>
53 {
55 
56 public:
57 
61 
62  typedef typename super::return_type return_type;
63  typedef T value_type;
65 
66  typedef typename super::nodes_type nodes_type;
67  typedef typename matrix_node<value_type>::type points_type;
68 
69  static const bool is_simplex = convex_type::is_simplex;
70  static const bool is_hypercube = convex_type::is_hypercube;
71 
72  static const uint32_type convexOrder = convex_type::nOrder;
73  static const uint32_type topological_dimension = convex_type::topological_dimension;
74 
75  typedef mpl::if_< mpl::bool_< is_simplex >,
77  Hypercube<Dim, Order, Dim> > conv_order_type;
78 
80 
81  static const uint32_type numPoints = conv_order_type::type::numPoints;
82 
83  static const uint32_type nbPtsPerVertex = conv_order_type::type::nbPtsPerVertex;
84  static const uint32_type nbPtsPerEdge = conv_order_type::type::nbPtsPerEdge;
85  static const uint32_type nbPtsPerFace = conv_order_type::type::nbPtsPerFace;
86  static const uint32_type nbPtsPerVolume = conv_order_type::type::nbPtsPerVolume;
87 
88  RefElem RefConv;
89 
90  typedef std::pair<uint16_type, uint16_type> range_type;
91  typedef std::vector< std::vector<size_type> > index_map_type;
92 
93 
95 
99 
101  :
102  super()
103  {
104  M_eid.resize( topological_dimension + 1 );
105  M_pt_to_entity.resize( numPoints );
106  }
108  :
109  super( np )
110  {
111  M_eid.resize( topological_dimension + 1 );
112  M_pt_to_entity.resize( numPoints );
113  }
115  :
116  super( psi ),
117  M_eid( psi.getEid() ),
118  M_pt_to_entity( psi.getPtE() )
119 
120  {}
121 
123 
124  ublas::matrix_range<nodes_type const> pointsByEntity( uint16_type e ) const
125  {
126  FEELPP_ASSERT( M_eid[e].size() )( e ).error( "no points defined on this entity" );
127 
128  return ublas::project( this->points(),
129  ublas::range( 0,Dim ),
130  ublas::range( *M_eid[e].begin(), *M_eid[e].rbegin() + 1 ) );
131  }
132 
133  range_type interiorRangeById( uint16_type e, uint16_type id ) const
134  {
135  uint16_type numEntities = 1;
136 
137  if ( e == 0 )
138  numEntities = convex_type::numVertices;
139 
140  else if ( e == 1 )
141  numEntities = convex_type::numEdges;
142 
143  else if ( e == 2 )
144  numEntities = convex_type::numFaces;
145 
146  uint16_type N = M_eid[e].size()/numEntities;
147 
148  return std::make_pair( *M_eid[e].begin() + id*N, *M_eid[e].begin() + ( id+1 )*N );
149  }
150 
151  ublas::matrix_range<nodes_type const> interiorPointsById( uint16_type e, uint16_type id ) const
152  {
153  range_type position = interiorRangeById( e, id );
154 
155  ublas::matrix_range<nodes_type const> G = ublas::project( this->points(),
156  ublas::range( 0, Dim ),
157  ublas::range( position.first, position.second ) );
158 
159  return G;
160  }
161 
162  uint32_type entityIds( int i, int j ) const
163  {
164  return M_eid[i][j];
165  }
166 
167  uint32_type numEntities( int i ) const
168  {
169  return M_eid[i].size();
170  }
171 
172  range_type const& pointToEntity( int p ) const
173  {
174  return M_pt_to_entity[p];
175  }
176 
177  //Returns the local indices of all the subentities that compose the entity
178  index_map_type entityToLocal ( uint16_type top_dim, uint16_type local_id, bool boundary = 0 )
179  {
180  index_map_type indices( top_dim+1 );
181 
182  if ( top_dim == 0 && boundary )
183  {
184  range_type pair = interiorRangeById( top_dim, local_id );
185 
186  indices[0].push_back( pair.first );
187  }
188 
189  else
190  {
191  if ( !boundary && top_dim > 0 )
192  {
193  if ( numEntities( top_dim ) != 0 )
194  indices[top_dim].push_back( local_id );
195  }
196 
197  else
198  {
199  if ( top_dim > 0 )
200  {
201  //For the time being, this definition works, but in a more general framework,
202  //a entity should be created here with the respective top_dim
203  //in order to retrieve the number of vertices, edges, faces.
204 
205  //number of vertices, number of edges and number of faces in the volume
206  std::map<uint16_type, uint16_type> numPointsInEntity;
207 
208  numPointsInEntity[0] = ( is_simplex )?( top_dim+1 ):( 2 + ( top_dim-1 )*( 2+3*( top_dim-2 ) ) );
209  numPointsInEntity[1] = ( 2 + ( top_dim-1 )*( 1+2*is_simplex + ( 5-4*is_simplex )*( top_dim-1 ) ) )/2;
210  numPointsInEntity[2] = 6-2*is_simplex;
211 
212  for ( uint16_type i=0; i < numPointsInEntity[0] ; i++ )
213  {
214  if ( top_dim == 1 )
215  indices[0].push_back( RefConv.e2p( local_id, i ) );
216 
217  else if ( top_dim == 2 )
218  indices[0].push_back( RefConv.f2p( local_id, i ) );
219  }
220 
221  if ( ( top_dim == 2 ) && ( M_eid[1].size() != 0 ) )
222  {
223  for ( uint16_type i=0; i < numPointsInEntity[1] ; i++ )
224  indices[1].push_back( RefConv.f2e( local_id, i ) );
225  }
226 
227  if ( top_dim == 3 )
228  {
229  for ( uint16_type k=0; k<numPointsInEntity.size(); k++ )
230  {
231  for ( uint16_type i=0; i < numPointsInEntity[k]; i++ )
232  {
233  indices[k].push_back( i );
234  }
235  }
236  }
237 
238  if ( M_eid[top_dim].size() )
239  {
240  indices[top_dim].push_back( local_id );
241  }
242  }
243  }
244  }
245 
246  return indices;
247  }
248 
249 
250  points_type pointsBySubEntity( uint16_type top_dim, uint16_type local_id, bool boundary = 0 )
251  {
252  index_map_type index_list = entityToLocal( top_dim, local_id, boundary );
253 
254  uint16_type matrix_size = 0;
255 
256  if ( index_list[0].size() != 0 )
257  matrix_size +=index_list[0].size()*nbPtsPerVertex;
258 
259  if ( ( top_dim>=1 ) && ( index_list[1].size() != 0 ) )
260  matrix_size +=index_list[1].size()*nbPtsPerEdge;
261 
262  if ( ( top_dim >=2 ) && ( index_list[2].size() != 0 ) )
263  matrix_size +=index_list[2].size()*nbPtsPerFace;
264 
265  if ( ( top_dim ==3 ) && ( index_list[3].size() != 0 ) )
266  matrix_size +=nbPtsPerVolume;
267 
268  points_type G ( Dim, matrix_size );
269 
270  for ( uint16_type i=0, p=0; i < top_dim+1; i++ )
271  {
272  if ( index_list[i].size() )
273  {
274  for ( uint16_type j=0; j < index_list[i].size(); j++ )
275  {
276  points_type aux = interiorPointsById( i, index_list[i][j] );
277 
278  ublas::subrange( G, 0, Dim, p, p+aux.size2() ) = aux;
279 
280  p+=aux.size2();
281  }
282  }
283  }
284 
285  return G;
286  }
287 
288  index_map_type getEid ()
289  {
290  return M_eid;
291  }
292 
293  std::vector<range_type> getPtE()
294  {
295  return M_pt_to_entity;
296  }
297 
298  void setEid ( index_map_type eid )
299  {
300  M_eid = eid;
301  }
302 
303  void setPtE ( std::vector<range_type> pt_ent )
304  {
305  M_pt_to_entity = pt_ent;
306  }
307 
308  void addToEid ( uint16_type p, uint16_type q )
309  {
310  M_eid[p].push_back( q );
311  }
312 
313  void addToPtE ( uint16_type p, range_type q )
314  {
315  M_pt_to_entity[p] = q;
316  }
317 
318  FEELPP_DEFINE_VISITABLE();
319 
320 private:
321 
322  index_map_type M_eid;
323  std::vector<range_type> M_pt_to_entity;
324 
325 };
326 }
327 #endif /* __PointSetInterpolation_H */

Generated on Sun Oct 20 2013 08:25:03 for Feel++ by doxygen 1.8.4