dune-fem 2.8.0
Loading...
Searching...
No Matches
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
15namespace 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 >
130 inline void LineSegmentSampler< GridPart >
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
155 DomainFieldType lower = std::numeric_limits< DomainFieldType >::min();
156 DomainFieldType upper = std::numeric_limits< DomainFieldType >::max();
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
182 lower = std::numeric_limits< DomainFieldType >::max();
183 upper = std::numeric_limits< DomainFieldType >::min();
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
Dune::Fem::Double abs(const Dune::Fem::Double &a)
Definition: double.hh:942
double max(const Dune::Fem::Double &v, const double p)
Definition: double.hh:965
double min(const Dune::Fem::Double &v, const double p)
Definition: double.hh:953
Definition: bindguard.hh:11
typename Impl::ConstLocalFunction< GridFunction >::Type ConstLocalFunction
Definition: const.hh:604
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
const GridPart & gridPart() const
obtain grid part on which the LineSegmentSampler works
Definition: linesegmentsampler.hh:98
GridPart GridPartType
Definition: linesegmentsampler.hh:40
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