dune-grid
2.2.0
|
00001 #ifndef DUNE_DGFS_HH 00002 #define DUNE_DGFS_HH 00003 00004 #include <dune/grid/common/intersection.hh> 00005 #include <dune/grid/sgrid.hh> 00006 #include "dgfparser.hh" 00007 00008 namespace Dune 00009 { 00010 // forward declaration 00011 // ------------------- 00012 00013 template< class GridImp, template< class > class IntersectionImp > 00014 class Intersection; 00015 00016 00017 template< int dim, int dimworld, class ctype > 00018 struct DGFGridFactory< SGrid< dim, dimworld, ctype > > 00019 { 00020 typedef SGrid< dim, dimworld, ctype > Grid; 00021 00022 const static int dimension = Grid::dimension; 00023 00024 typedef MPIHelper::MPICommunicator MPICommunicatorType; 00025 00026 private: 00027 typedef FieldVector< double, dimension > Point; 00028 00029 typedef dgf::BoundaryDomBlock BoundaryDomainBlock; 00030 00031 public: 00032 explicit DGFGridFactory ( std::istream &input, 00033 MPICommunicatorType comm = MPIHelper::getCommunicator() ) 00034 { 00035 generate( input, comm ); 00036 } 00037 00038 explicit DGFGridFactory ( const std::string &filename, 00039 MPICommunicatorType comm = MPIHelper::getCommunicator() ) 00040 { 00041 std::ifstream input( filename.c_str() ); 00042 generate( input, comm ); 00043 } 00044 00045 Grid *grid() const 00046 { 00047 return grid_; 00048 } 00049 00050 template< class Intersection > 00051 bool wasInserted ( const Intersection &intersection ) const 00052 { 00053 return false; 00054 } 00055 00056 template< class Intersection > 00057 int boundaryId ( const Intersection &intersection ) const 00058 { 00059 if( boundaryDomainBlock_->isactive() ) 00060 { 00061 std::vector< Point > corners; 00062 getCorners( intersection.geometry(), corners ); 00063 const dgf::DomainData *data = boundaryDomainBlock_->contains( corners ); 00064 if( data ) 00065 return data->id(); 00066 else 00067 return intersection.indexInInside(); 00068 } 00069 else 00070 return intersection.indexInInside(); 00071 } 00072 00073 template< int codim > 00074 int numParameters () const 00075 { 00076 return 0; 00077 } 00078 00079 // return true if boundary parameters found 00080 bool haveBoundaryParameters () const 00081 { 00082 return boundaryDomainBlock_->isactive(); 00083 } 00084 00085 template < class GG, template < class > class II > 00086 const typename DGFBoundaryParameter::type & 00087 boundaryParameter ( const Intersection< GG, II > & intersection ) const 00088 { 00089 if( haveBoundaryParameters() ) 00090 { 00091 std::vector< Point > corners; 00092 getCorners( intersection.geometry(), corners ); 00093 const dgf::DomainData *data = boundaryDomainBlock_->contains( corners ); 00094 if( data ) 00095 return data->parameter(); 00096 else 00097 return DGFBoundaryParameter::defaultValue(); 00098 } 00099 else 00100 return DGFBoundaryParameter::defaultValue(); 00101 } 00102 00103 00104 template< class Entity > 00105 std::vector< double > ¶meter ( const Entity &entity ) 00106 { 00107 return emptyParam; 00108 } 00109 00110 private: 00111 void generate( std::istream &gridin, MPICommunicatorType comm ); 00112 00113 template< class Geometry > 00114 static void getCorners ( const Geometry &geometry, std::vector< Point > &corners ) 00115 { 00116 corners.resize( geometry.corners() ); 00117 for( int i = 0; i < geometry.corners(); ++i ) 00118 { 00119 const typename Geometry::GlobalCoordinate corner = geometry.corner( i ); 00120 for( int j = 0; j < dimension; ++j ) 00121 corners[ i ][ j ] = corner[ j ]; 00122 } 00123 } 00124 00125 Grid *grid_; 00126 dgf::BoundaryDomBlock *boundaryDomainBlock_; 00127 std::vector< double > emptyParam; 00128 }; 00129 00130 00131 00132 template< int dim, int dimworld, class ctype > 00133 inline void DGFGridFactory< SGrid< dim, dimworld, ctype > > 00134 ::generate ( std::istream &gridin, MPICommunicatorType comm ) 00135 { 00136 dgf::IntervalBlock intervalBlock( gridin ); 00137 00138 if( !intervalBlock.isactive() ) 00139 DUNE_THROW( DGFException, "SGrid can only be created from an interval block." ); 00140 00141 if( intervalBlock.numIntervals() != 1 ) 00142 DUNE_THROW( DGFException, "SGrid can only handle 1 interval block." ); 00143 00144 if( intervalBlock.dimw() != dim ) 00145 { 00146 DUNE_THROW( DGFException, 00147 "Cannot read an interval of dimension " << intervalBlock.dimw() 00148 << "into a SGrid< " << dim << ", " << dimworld << " >." ); 00149 } 00150 00151 const dgf::IntervalBlock::Interval &interval = intervalBlock.get( 0 ); 00152 00153 FieldVector< double, dimension > lower, upper; 00154 FieldVector< int, dimension > anz; 00155 for( int i = 0; i < dimension; ++i ) 00156 { 00157 lower[ i ] = interval.p[ 0 ][ i ]; 00158 upper[ i ] = interval.p[ 1 ][ i ]; 00159 anz[ i ] = interval.n[ i ]; 00160 } 00161 00162 grid_ = new Grid( anz, lower, upper ); 00163 00164 boundaryDomainBlock_ = new dgf::BoundaryDomBlock( gridin, dimension ); 00165 } 00166 00167 00168 00169 template< int dim, int dimworld, class ctype > 00170 struct DGFGridInfo< SGrid< dim, dimworld, ctype > > 00171 { 00172 static int refineStepsForHalf () 00173 { 00174 return 1; 00175 } 00176 00177 static double refineWeight () 00178 { 00179 return 1.0 / double( 1 << dim ); 00180 } 00181 }; 00182 00183 } 00184 #endif // #ifndef DUNE_DGFS_HH