1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_GRID_FILE_IO_DGFPARSER_DGFONED_HH
4 #define DUNE_GRID_FILE_IO_DGFPARSER_DGFONED_HH
5 
6 //- C++ includes
7 #include <algorithm>
8 #include <fstream>
9 #include <iostream>
10 #include <istream>
11 #include <vector>
12 
13 //- dune-common includes
14 #include <dune/common/exceptions.hh>
15 
16 //- dune-grid includes
17 #include <dune/grid/common/intersection.hh>
18 #include <dune/grid/onedgrid.hh>
19 
20 //- local includes
21 #include "dgfparser.hh"
22 
23 
24 namespace
25 {
26   // helper method used below
getfirst(std::vector<double> v)27   double getfirst ( std::vector< double > v )
28   {
29     return v[ 0 ];
30   }
31 }  // end anonymous namespace
32 
33 
34 
35 namespace Dune
36 {
37 
38   // DGFGridInfo
39   // -----------
40 
41   template< >
42   struct DGFGridInfo< OneDGrid >
43   {
refineStepsForHalfDune::DGFGridInfo44     static int refineStepsForHalf ()
45     {
46       return 1;
47     }
48 
refineWeightDune::DGFGridInfo49     static double refineWeight ()
50     {
51       return 0.5;
52     }
53   };
54 
55 
56 
57   // DGFGridFactory< OneDGrid >
58   // --------------------------
59 
60   template< >
61   struct DGFGridFactory< OneDGrid >
62   {
63     /** \brief grid type */
64     typedef OneDGrid Grid;
65     /** \brief grid dimension */
66     const static int dimension = Grid::dimension;
67     /** \brief MPI communicator type */
68     typedef MPIHelper::MPICommunicator MPICommunicatorType;
69 
70     /** \brief constructor taking istream */
DGFGridFactoryDune::DGFGridFactory71     explicit DGFGridFactory ( std::istream &input,
72                               MPICommunicatorType comm = MPIHelper::getCommunicator() )
73       : grid_( 0 ),
74         emptyParameters_( 0 )
75     {
76       generate( input, comm );
77     }
78 
79     /** \brief constructor taking filename */
DGFGridFactoryDune::DGFGridFactory80     explicit DGFGridFactory ( const std::string &filename,
81                               MPICommunicatorType comm = MPIHelper::getCommunicator() )
82       : grid_( 0 ),
83         emptyParameters_( 0 )
84     {
85       std::ifstream input( filename.c_str() );
86       generate( input, comm );
87     }
88 
89     /** \brief get grid */
gridDune::DGFGridFactory90     Grid *grid () const
91     {
92       return grid_;
93     }
94 
95     /** \brief always returns false */
96     template< class GG, class II >
wasInsertedDune::DGFGridFactory97     bool wasInserted ( const Dune::Intersection< GG, II > &intersection ) const
98     {
99       return false;
100     }
101 
102     template< class GG, class II >
boundaryIdDune::DGFGridFactory103     int boundaryId ( const Dune::Intersection< GG, II > &intersection ) const
104     {
105       // OneDGrid returns boundary segment index;
106       // we return the index as the method boundaryId is deprecated
107       return intersection.boundarySegmentIndex();
108     }
109 
110     /** \brief OneDGrid does not support parameters, returns 0 */
111     template< class Entity >
numParametersDune::DGFGridFactory112     int numParameters ( const Entity & ) const
113     {
114       return 0;
115     }
116 
117     /** \brief OneDGrid does not support parameters, returns 0 */
118     template< int codim >
numParametersDune::DGFGridFactory119     int numParameters () const
120     {
121       return 0;
122     }
123 
124     template< class Entity >
parameterDune::DGFGridFactory125     std::vector< double >& parameter ( const Entity &entity )
126     {
127       return parameter< Entity::codimension >( entity );
128     }
129 
130     /** \brief return empty vector */
131     template< int codim >
parameterDune::DGFGridFactory132     std::vector< double > &parameter ( [[maybe_unused]] const typename Grid::Codim< codim >::Entity &element )
133     {
134       return emptyParameters_;
135     }
136 
137     /** \brief OneDGrid does not support boundary parameters */
haveBoundaryParametersDune::DGFGridFactory138     bool haveBoundaryParameters () const
139     {
140       return false;
141     }
142 
143     /** \brief return invalid default value */
144     template< class GG, class II >
boundaryParameterDune::DGFGridFactory145     const DGFBoundaryParameter::type &boundaryParameter ( [[maybe_unused]] const Dune::Intersection< GG, II > &intersection ) const
146     {
147       return DGFBoundaryParameter::defaultValue();
148     }
149 
150   private:
151     // generate grid
152     void generate ( std::istream &input, MPICommunicatorType comm );
153 
154     Grid *grid_;
155     std::vector< double > emptyParameters_;
156   };
157 
158 
159 
160   // Implementation of DGFGridFactory< OneDGrid >
161   // --------------------------------------------
162 
generate(std::istream & input,MPICommunicatorType comm)163   inline void DGFGridFactory< OneDGrid >::generate ( std::istream &input, [[maybe_unused]] MPICommunicatorType comm )
164   {
165     // try to open interval block
166     dgf::IntervalBlock intervalBlock( input );
167 
168     // try to open vertex block
169     int dimensionworld = Grid::dimensionworld;
170     dgf::VertexBlock vertexBlock( input, dimensionworld );
171 
172     // check at least one block is active
173     if( !( vertexBlock.isactive() || intervalBlock.isactive() ))
174       DUNE_THROW( DGFException, "No readable block found" );
175 
176     std::vector< std::vector< double > > vertices;
177 
178     // read vertices first
179     if( vertexBlock.isactive() )
180     {
181       int nparameter = 0;
182       std::vector< std::vector< double > > parameter;
183       vertexBlock.get( vertices, parameter, nparameter );
184 
185       if( nparameter > 0 )
186         std::cerr << "Warning: vertex parameters will be ignored" << std::endl;
187     }
188 
189     // get vertices from interval block
190     if ( intervalBlock.isactive() )
191     {
192       if( intervalBlock.dimw() != dimensionworld )
193       {
194         DUNE_THROW( DGFException, "Error: wrong coordinate dimension in interval block \
195                                    (got "                                                                                         << intervalBlock.dimw() << ", expected " << dimensionworld << ")" );
196       }
197 
198       int nintervals = intervalBlock.numIntervals();
199       for( int i = 0; i < nintervals; ++i )
200         intervalBlock.getVtx( i, vertices );
201     }
202 
203     // copy to vector of doubles
204     std::vector< double > vtx( vertices.size() );
205     transform( vertices.begin(), vertices.end(), vtx.begin(), getfirst );
206 
207     // remove duplicates
208     std::sort( vtx.begin(), vtx.end() );
209     std::vector< double >::iterator it = std::unique( vtx.begin(), vtx.end() );
210     vtx.erase( it, vtx.end() );
211     if( vertices.size() != vtx.size() )
212       std::cerr << "Warning: removed duplicate vertices" << std::endl;
213 
214     // create grid
215     grid_ = new OneDGrid( vtx );
216   }
217 
218 } // end namespace Dune
219 
220 #endif // #ifndef DUNE_GRID_FILE_IO_DGFPARSER_DGFONED_HH
221