dune-fem  2.8-git
linesegmentsampler.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_LINESEGMENTSAMPLER_HH
2 #define DUNE_FEM_LINESEGMENTSAMPLER_HH
3 
4 #include <limits>
5 #include <vector>
6 #include <cmath>
7 #include <type_traits>
8 
9 #include <dune/common/fvector.hh>
10 
11 #include <dune/geometry/referenceelements.hh>
12 
14 
15 namespace Dune
16 {
17 
18  namespace Fem
19  {
20 
21  // LineSegmentSampler
22  // ------------------
23 
37  template< class GridPart >
39  {
40  typedef GridPart GridPartType;
41 
42  typedef typename GridPartType::GridType::ctype DomainFieldType;
43 
44  static const int dimDomain = GridPartType::dimensionworld;
45  static const int dimGrid = GridPartType::dimension;
46 
47  typedef FieldVector< DomainFieldType, dimDomain > DomainType;
48  typedef FieldVector< DomainFieldType, dimGrid > LocalDomainType;
49 
50  private:
51  static_assert( dimDomain == dimGrid, "LineSegmentSampler supports only flat grids." );
52 
53  template< class Vector >
54  struct Reduce;
55 
56  typedef Dune::ReferenceElements< DomainFieldType, dimGrid > ReferenceElements;
57 
58  public:
69  LineSegmentSampler ( const GridPart &gridPart, const DomainType &left, const DomainType &right )
70  : gridPart_( gridPart ), left_( left ), right_( right )
71  {}
72 
83  template< class GridFunction >
84  void operator() ( const GridFunction &f, std::vector< typename GridFunction::RangeType > &samples ) const;
85 
95  void samplePoints( std::vector< DomainType > &points ) const;
96 
98  const GridPart &gridPart () const { return gridPart_; }
99 
100  private:
101  const GridPart &gridPart_;
102  DomainType left_, right_;
103  };
104 
105 
106 
107  // LineSegmentSampler::Reduce
108  // --------------------------
109 
110  template< class GridPart >
111  template< class Vector >
112  struct LineSegmentSampler< GridPart >::Reduce
113  {
114  Vector operator() ( const Vector &a, const Vector &b ) const
115  {
116  Vector c;
117  for( int k = 0; k < Vector::dimension; ++k )
118  c[ k ] = std::min( a[ k ], b[ k ] );
119  return c;
120  }
121  };
122 
123 
124 
125  // Implementation of LineSegmentSampler
126  // ------------------------------------
127 
128  template< class GridPart >
129  template< class GridFunction >
131  ::operator() ( const GridFunction &f, std::vector< typename GridFunction::RangeType > &samples ) const
132  {
133  typedef typename GridPartType::template Codim< 0 >::IteratorType IteratorType;
134  typedef typename GridPartType::template Codim< 0 >::EntityType EntityType;
135 
136  const int numSamples = samples.size();
137  if( numSamples < 2 )
138  DUNE_THROW( InvalidStateException, "LineSegmentSampler cannot sample less than 2 points." );
139  DomainType ds = right_ - left_;
140  ds /= DomainFieldType( numSamples - 1 );
141 
142  typedef typename GridFunction::RangeFieldType RangeFieldType;
143  const RangeFieldType invalid
144  = std::numeric_limits< RangeFieldType >::quiet_NaN();
145  for( int i = 0; i < numSamples; ++i )
146  samples[ i ] = typename GridFunction::RangeType( invalid );
147 
149  const IteratorType end = gridPart().template end< 0 >();
150  for( IteratorType it = gridPart().template begin< 0 >(); it != end; ++it )
151  {
152  const EntityType &entity = *it;
153  const typename EntityType::Geometry &geometry = entity.geometry();
154 
157 
158  auto &refElement = ReferenceElements::general( geometry.type() );
159  const int numFaces = refElement.size( 1 );
160  for( int face = 0; face < numFaces; ++face )
161  {
162  const LocalDomainType &center = refElement.position( face, 1 );
163 
164  DomainType normal;
165  const LocalDomainType &refNormal = refElement.integrationOuterNormal( face );
166  geometry.jacobianInverseTransposed( center ).mv( refNormal, normal );
167 
168  const DomainFieldType nds = normal * ds;
169  const DomainFieldType ncl = normal * (geometry.global( center ) - left_);
170  if( std::abs( nds ) > 1e-8 )
171  {
172  // ds is not parallel to the face
173  const DomainFieldType alpha = ncl / nds;
174  if( nds > 0 )
175  upper = std::min( upper, alpha );
176  else
177  lower = std::max( lower, alpha );
178  }
179  else if( ncl < -1e-8 )
180  {
181  // ds is parallel to the face and the line lies on the outside
184  }
185  }
186 
187  if( lower <= upper )
188  {
189  lf.bind( entity );
190  const int i_upper = std::min( int( std::floor( upper + 1e-8 ) ), numSamples - 1 );
191  const int i_lower = std::max( int( std::ceil( lower - 1e-8 ) ), 0 );
192  for( int i = i_lower; i <= i_upper; ++i )
193  {
194  DomainType xi = left_;
195  xi.axpy( DomainFieldType( i ), ds );
196  lf.evaluate( geometry.local( xi ), samples[ i ] );
197  }
198  lf.unbind();
199  }
200  }
201 
202  typedef Reduce< typename GridFunction::RangeType > Op;
203  gridPart().comm().template allreduce< Op >( &(samples[ 0 ]), numSamples );
204 
205  bool valid = true;
206 
207  // only use isnan when field type is a floating point
208  if constexpr ( std::is_floating_point< RangeFieldType >::value )
209  {
210  for( int i = 0; i < numSamples; ++i )
211  {
212  for( int d=0; d<GridFunction::RangeType::dimension; ++d )
213  {
214  valid &= ! std::isnan( samples[ i ][ d ] );
215  }
216  }
217  }
218 
219  if( !valid )
220  DUNE_THROW( InvalidStateException, "LineSegmentSampler could not find all samples." );
221  }
222 
223 
224  template< class GridPart >
226  :: samplePoints ( std::vector< DomainType > &points ) const
227  {
228  const int numSamples = points.size();
229  if( numSamples < 2 )
230  DUNE_THROW( InvalidStateException, "LineSegmentSampler cannot sample less than 2 points." );
231  DomainType ds = right_ - left_;
232  ds /= DomainFieldType( numSamples - 1 );
233 
234  for( int i = 0; i < numSamples; ++i )
235  {
236  points[ i ] = left_;
237  points[ i ].axpy( DomainFieldType( i ), ds );
238  }
239  }
240 
241  } // namespace Fem
242 
243 } // namespace Dune
244 
245 #endif // #ifndef DUNE_FEM_LINESEGMENTSAMPLER_HH
Definition: bindguard.hh:11
Double abs(const Double &a)
Definition: double.hh:871
typename Impl::ConstLocalFunction< GridFunction >::Type ConstLocalFunction
Definition: const.hh:517
static constexpr T max(T a)
Definition: utility.hh:77
static constexpr T min(T a)
Definition: utility.hh:93
samples values of a discrete function along a given line segment
Definition: linesegmentsampler.hh:39
void samplePoints(std::vector< DomainType > &points) const
returns sampling points
Definition: linesegmentsampler.hh:226
GridPart GridPartType
Definition: linesegmentsampler.hh:40
const GridPart & gridPart() const
obtain grid part on which the LineSegmentSampler works
Definition: linesegmentsampler.hh:98
static const int dimDomain
Definition: linesegmentsampler.hh:44
void operator()(const GridFunction &f, std::vector< typename GridFunction::RangeType > &samples) const
sample a given function
Definition: linesegmentsampler.hh:131
LineSegmentSampler(const GridPart &gridPart, const DomainType &left, const DomainType &right)
constructor
Definition: linesegmentsampler.hh:69
static const int dimGrid
Definition: linesegmentsampler.hh:45
FieldVector< DomainFieldType, dimDomain > DomainType
Definition: linesegmentsampler.hh:47
GridPartType::GridType::ctype DomainFieldType
Definition: linesegmentsampler.hh:42
FieldVector< DomainFieldType, dimGrid > LocalDomainType
Definition: linesegmentsampler.hh:48