29 #ifndef __StructuredGrid_H
30 #define __StructuredGrid_H 1
35 #include <boost/multi_array.hpp>
36 #include <boost/array.hpp>
60 min = ublas::scalar_vector<double>( Dim, -1 );
61 max = ublas::scalar_vector<double>( Dim, 1 );
69 bbox( ublas::vector<node_type>
const& __cp )
73 eps( ublas::scalar_vector<double>( min.size(), 1e-10 ) )
78 for (
int __i = 1; __i < __cp.size(); ++__i )
80 for (
int __d = 0; __d < Dim; ++__d )
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];
97 bbox( std::vector<boost::shared_ptr<bbox<Dim> > >
const& bboxes )
101 eps( ublas::scalar_vector<double>( min.size(), 1e-10 ) )
103 min = bboxes[0]->min;
104 max = bboxes[0]->max;
106 for (
int __i = 1; __i < bboxes.size(); ++__i )
108 for (
int __d = 0; __d < Dim; ++__d )
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];
125 void enlarge(
double e )
127 min -= ublas::scalar_vector<double>( min.size(), e );
128 max += ublas::scalar_vector<double>( min.size(), e );
142 template<u
int16_type Dim, u
int16_type Order = 1>
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;
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;
177 StructuredGrid( bbox<Dim>
const& __bbox, shape_type
const& __shape )
183 std::for_each( M_coord.begin(), M_coord.end(), lambda::bind( &node_type::resize,
188 node_type __pt( Dim );
190 for ( index_type i = 0; i < M_shape[0]; ++i )
193 __pt[0] = __bbox.min[0]+i*( __bbox.max[0]-__bbox.min[0] )/( M_shape[0]-1 );
195 for ( index_type j = 0; j < M_shape[1]; ++j )
198 __pt[1] = __bbox.min[1]+j*( __bbox.max[1]-__bbox.min[1] )/( M_shape[1]-1 );
200 for ( index_type k = 0; k < M_shape[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;
232 return M_coord.num_elements();
241 shape_type
const& shape()
const
253 uint extent( uint __i )
const
264 node_type& operator()( idx_type
const& __idx )
266 return M_coord( __idx );
275 node_type operator()( idx_type
const& __idx )
const
277 return M_coord( __idx );
298 void save( std::string
const& prefix,
299 std::vector<boost::shared_ptr<f_type> >
const& __f );
315 template<u
int16_type Dim, u
int16_type Order>
318 std::vector<boost::shared_ptr<f_type> >
const& __f )
320 std::ofstream __ofs( ( prefix +
".geo" ).c_str() );
327 << std::setw( 8 ) << 0 <<
"\n";
333 << std::setw( 8 ) << M_shape[0]
334 << std::setw( 8 ) << M_shape[1]
335 << std::setw( 8 ) << M_shape[2];
340 for ( index_type i = 0; i < M_shape[0]; ++i )
344 for ( index_type j = 0; j < M_shape[1]; ++j )
348 for ( index_type k = 0; k < M_shape[2]; ++k )
352 if ( __count ++ % 6 == 0 )
355 __ofs.precision( 5 );
356 __ofs.setf( std::ios::scientific );
357 __ofs << std::setw( 12 ) << M_coord( __idx )[0];
366 for ( index_type i = 0; i < M_shape[0]; ++i )
370 for ( index_type j = 0; j < M_shape[1]; ++j )
374 for ( index_type k = 0; k < M_shape[2]; ++k )
378 if ( __count ++ % 6 == 0 )
381 __ofs.precision( 5 );
382 __ofs.setf( std::ios::scientific );
383 __ofs << std::setw( 12 ) << M_coord( __idx )[1];
392 for ( index_type i = 0; i < M_shape[0]; ++i )
396 for ( index_type j = 0; j < M_shape[1]; ++j )
400 for ( index_type k = 0; k < M_shape[2]; ++k )
404 if ( __count ++ % 6 == 0 )
407 __ofs.precision( 5 );
408 __ofs.setf( std::ios::scientific );
409 __ofs << std::setw( 12 ) << M_coord( __idx )[2];
417 __ofs.open( ( prefix +
".case" ).c_str() );
425 for (
int i = 0; i < __f.size(); ++i )
427 __ofs <<
"scalar per node: f" << i <<
" f" << i <<
".001\n";
432 for (
int l = 0; l < __f.size(); ++l )
434 std::ostringstream __ostr;
435 __ostr <<
"f" << l <<
".001";
436 __ofs.open( __ostr.str().c_str() );
438 __ofs.precision( 5 );
439 __ofs.setf( std::ios::scientific );
440 __ofs <<
"function\n"
445 for ( index_type i = 0; i < M_shape[0]; ++i )
449 for ( index_type j = 0; j < M_shape[1]; ++j )
453 for ( index_type k = 0; k < M_shape[2]; ++k )
457 if ( __count ++ % 6 == 0 )
460 __ofs << std::setw( 12 ) << __f[l]->operator()( __idx );