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
structuredgrid.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 Imaging project.
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2005-02-28
7 
8  Copyright (C) 2005,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 __StructuredGrid_H
30 #define __StructuredGrid_H 1
31 #include <iostream>
32 #include <fstream>
33 #include <iomanip>
34 
35 #include <boost/multi_array.hpp>
36 #include <boost/array.hpp>
37 #include <feel/feelalg/glas.hpp>
38 
39 namespace Feel
40 {
42 
49 template<uint Dim>
50 struct bbox
51 {
55  bbox()
56  :
57  min( Dim ),
58  max( Dim )
59  {
60  min = ublas::scalar_vector<double>( Dim, -1 );
61  max = ublas::scalar_vector<double>( Dim, 1 );
62  }
63 
69  bbox( ublas::vector<node_type> const& __cp )
70  :
71  min( Dim ),
72  max( Dim ),
73  eps( ublas::scalar_vector<double>( min.size(), 1e-10 ) )
74  {
75  min = __cp( 0 );
76  max = __cp( 0 );
77 
78  for ( int __i = 1; __i < __cp.size(); ++__i )
79  {
80  for ( int __d = 0; __d < Dim; ++__d )
81  {
82  min[__d] = ( __cp( __i )[__d] > min[__d] )?min[__d]: __cp( __i )[__d];
83  max[__d] = ( __cp( __i )[__d] < max[__d] )?max[__d]: __cp( __i )[__d];
84  }
85  }
86 
87  // enlarge a tad the bounding box
88  min -= eps;
89  max += eps;
90  }
91 
97  bbox( std::vector<boost::shared_ptr<bbox<Dim> > > const& bboxes )
98  :
99  min( Dim ),
100  max( Dim ),
101  eps( ublas::scalar_vector<double>( min.size(), 1e-10 ) )
102  {
103  min = bboxes[0]->min;
104  max = bboxes[0]->max;
105 
106  for ( int __i = 1; __i < bboxes.size(); ++__i )
107  {
108  for ( int __d = 0; __d < Dim; ++__d )
109  {
110  min[__d] = ( bboxes[__i]->min[__d] > min[__d] )?min[__d]: bboxes[ __i ]->min[__d];
111  max[__d] = ( bboxes[__i]->max[__d] < max[__d] )?max[__d]: bboxes[ __i ]->max[__d];
112  }
113  }
114 
115  // enlarge a tad the bounding box
116  min -= eps;
117  max += eps;
118  }
119 
125  void enlarge( double e )
126  {
127  min -= ublas::scalar_vector<double>( min.size(), e );
128  max += ublas::scalar_vector<double>( min.size(), e );
129  }
130  node_type min;
131  node_type max;
132  node_type eps;
133 };
135 
142 template<uint16_type Dim, uint16_type Order = 1>
144 {
145 public:
146 
147 
151 
152  typedef boost::multi_array<node_type, Dim> coord_type;
153  typedef boost::multi_array<double, Dim> f_type;
154  typedef typename coord_type::index index_type;
155  typedef boost::array<index_type, Dim> shape_type;
156  typedef boost::array<index_type, Dim> idx_type;
157 
158  typedef typename mpl::if_<mpl::equal<mpl::int_<Dim>, mpl::int_<1> >,
159  mpl::identity<GeoElement1D<Dim, LinearLine, DefMarkerCommon> >,
160  mpl::if_<mpl::equal<mpl::int_<Dim>, mpl::int_<2> >,
161  mpl::identity<GeoElement2D<Dim, LinearQuad, DefMarkerCommon> >,
162  mpl::if_<mpl::equal<mpl::int_<Dim>, mpl::int_<3> >,
163  mpl::identity<GeoElement3D<Dim, LinearHexa, DefMarkerCommon> > >::type::type element_type;
164  typedef typename mpl::if_<mpl::equal<mpl::int_<Dim>, mpl::int_<1> >,
165  mpl::identity<GeoElement0D<Dim, DefMarkerCommon> >,
166  mpl::if_<mpl::equal<mpl::int_<Dim>, mpl::int_<2> >,
167  mpl::identity<GeoElement1D<Dim, LinearLine, DefMarkerCommon> >,
168  mpl::if_<mpl::equal<mpl::int_<Dim>, mpl::int_<3> >,
169  mpl::identity<GeoElement2D<Dim, LinearQuad, DefMarkerCommon> > >::type::type face_type;
170 
172 
176 
177  StructuredGrid( bbox<Dim> const& __bbox, shape_type const& __shape )
178  :
179  M_shape( __shape ),
180  M_coord( __shape )
181  {
182 #if 0
183  std::for_each( M_coord.begin(), M_coord.end(), lambda::bind( &node_type::resize,
184  lambda::_1,
185  Dim ) );
186 #else
187  idx_type __idx;
188  node_type __pt( Dim );
189 
190  for ( index_type i = 0; i < M_shape[0]; ++i )
191  {
192  __idx[0] = i;
193  __pt[0] = __bbox.min[0]+i*( __bbox.max[0]-__bbox.min[0] )/( M_shape[0]-1 );
194 
195  for ( index_type j = 0; j < M_shape[1]; ++j )
196  {
197  __idx[1] = j;
198  __pt[1] = __bbox.min[1]+j*( __bbox.max[1]-__bbox.min[1] )/( M_shape[1]-1 );
199 
200  for ( index_type k = 0; k < M_shape[2]; ++k )
201  {
202  __idx[2] = k;
203  __pt[2] = __bbox.min[2]+k*( __bbox.max[2]-__bbox.min[2] )/( M_shape[2]-1 );
204  M_coord( __idx ).resize( Dim );
205  M_coord( __idx ) = __pt;
206  }
207  }
208  }
209 
210 #endif
211  }
212  ~StructuredGrid()
213  {}
214 
216 
220 
221 
223 
227 
230  size_type numElements() const
231  {
232  return M_coord.num_elements();
233  }
234 
241  shape_type const& shape() const
242  {
243  return M_shape;
244  }
245 
253  uint extent( uint __i ) const
254  {
255  return M_shape[__i];
256  }
257 
264  node_type& operator()( idx_type const& __idx )
265  {
266  return M_coord( __idx );
267  }
268 
275  node_type operator()( idx_type const& __idx ) const
276  {
277  return M_coord( __idx );
278  }
279 
281 
285 
286 
288 
292 
298  void save( std::string const& prefix,
299  std::vector<boost::shared_ptr<f_type> > const& __f );
300 
302 
303 protected:
304 
305  shape_type M_shape;
306 
307  coord_type M_coord;
308 
309 private:
310 
311  StructuredGrid();
312  StructuredGrid( StructuredGrid const & );
313 
314 };
315 template<uint16_type Dim, uint16_type Order>
316 void
317 StructuredGrid<Dim, Order>::save( std::string const& prefix,
318  std::vector<boost::shared_ptr<f_type> > const& __f )
319 {
320  std::ofstream __ofs( ( prefix + ".geo" ).c_str() );
321  __ofs <<
322  "Structured grid\n"
323  "Structured grid\n"
324  "node id off\n"
325  "element id off\n"
326  "coordinates\n"
327  << std::setw( 8 ) << 0 << "\n";
328 
329  __ofs <<
330  "part 1\n"
331  "structured block\n"
332  "block\n"
333  << std::setw( 8 ) << M_shape[0]
334  << std::setw( 8 ) << M_shape[1]
335  << std::setw( 8 ) << M_shape[2];
336 
337  idx_type __idx;
338  size_t __count = 0;
339 
340  for ( index_type i = 0; i < M_shape[0]; ++i )
341  {
342  __idx[0] = i;
343 
344  for ( index_type j = 0; j < M_shape[1]; ++j )
345  {
346  __idx[1] = j;
347 
348  for ( index_type k = 0; k < M_shape[2]; ++k )
349  {
350  __idx[2] = k;
351 
352  if ( __count ++ % 6 == 0 )
353  __ofs << "\n";
354 
355  __ofs.precision( 5 );
356  __ofs.setf( std::ios::scientific );
357  __ofs << std::setw( 12 ) << M_coord( __idx )[0];
358  }
359  }
360  }
361 
362  //if ( __count-1 % 6 != 0 )
363  // __ofs << "\n";
364  __count = 0;
365 
366  for ( index_type i = 0; i < M_shape[0]; ++i )
367  {
368  __idx[0] = i;
369 
370  for ( index_type j = 0; j < M_shape[1]; ++j )
371  {
372  __idx[1] = j;
373 
374  for ( index_type k = 0; k < M_shape[2]; ++k )
375  {
376  __idx[2] = k;
377 
378  if ( __count ++ % 6 == 0 )
379  __ofs << "\n";
380 
381  __ofs.precision( 5 );
382  __ofs.setf( std::ios::scientific );
383  __ofs << std::setw( 12 ) << M_coord( __idx )[1];
384  }
385  }
386  }
387 
388  //if ( __count-1 % 6 != 0 )
389  // __ofs << "\n";
390  __count = 0;
391 
392  for ( index_type i = 0; i < M_shape[0]; ++i )
393  {
394  __idx[0] = i;
395 
396  for ( index_type j = 0; j < M_shape[1]; ++j )
397  {
398  __idx[1] = j;
399 
400  for ( index_type k = 0; k < M_shape[2]; ++k )
401  {
402  __idx[2] = k;
403 
404  if ( __count ++ % 6 == 0 )
405  __ofs << "\n";
406 
407  __ofs.precision( 5 );
408  __ofs.setf( std::ios::scientific );
409  __ofs << std::setw( 12 ) << M_coord( __idx )[2];
410  }
411  }
412  }
413 
414  __ofs << "\n";
415 
416  __ofs.close();
417  __ofs.open( ( prefix + ".case" ).c_str() );
418  __ofs <<
419  "FORMAT\n"
420  "type: ensight\n"
421  "GEOMETRY\n"
422  "model: toto.geo\n"
423  "VARIABLE\n";
424 
425  for ( int i = 0; i < __f.size(); ++i )
426  {
427  __ofs << "scalar per node: f" << i << " f" << i << ".001\n";
428  }
429 
430  __ofs.close();
431 
432  for ( int l = 0; l < __f.size(); ++l )
433  {
434  std::ostringstream __ostr;
435  __ostr << "f" << l << ".001";
436  __ofs.open( __ostr.str().c_str() );
437 
438  __ofs.precision( 5 );
439  __ofs.setf( std::ios::scientific );
440  __ofs << "function\n"
441  << "part 1\n"
442  << "block";
443  __count = 0;
444 
445  for ( index_type i = 0; i < M_shape[0]; ++i )
446  {
447  __idx[0] = i;
448 
449  for ( index_type j = 0; j < M_shape[1]; ++j )
450  {
451  __idx[1] = j;
452 
453  for ( index_type k = 0; k < M_shape[2]; ++k )
454  {
455  __idx[2] = k;
456 
457  if ( __count ++ % 6 == 0 )
458  __ofs << "\n";
459 
460  __ofs << std::setw( 12 ) << __f[l]->operator()( __idx );
461  }
462  }
463  }
464 
465  __ofs.close();
466  }
467 } // end StructuredGrid<>::save
468 
469 }
470 #endif /* __StructuredGrid_H */

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