1#ifndef DUNE_FEM_LINESEGMENTSAMPLER_HH
2#define DUNE_FEM_LINESEGMENTSAMPLER_HH
9#include <dune/common/fvector.hh>
11#include <dune/geometry/referenceelements.hh>
37 template<
class Gr
idPart >
44 static const int dimDomain = GridPartType::dimensionworld;
45 static const int dimGrid = GridPartType::dimension;
47 typedef FieldVector< DomainFieldType, dimDomain >
DomainType;
51 static_assert(
dimDomain ==
dimGrid,
"LineSegmentSampler supports only flat grids." );
53 template<
class Vector >
56 typedef Dune::ReferenceElements< DomainFieldType, dimGrid > ReferenceElements;
70 : gridPart_(
gridPart ), left_( left ), right_( right )
83 template<
class Gr
idFunction >
84 void operator() (
const GridFunction &f, std::vector< typename GridFunction::RangeType > &samples )
const;
95 void samplePoints( std::vector< DomainType > &points )
const;
98 const GridPart &
gridPart ()
const {
return gridPart_; }
101 const GridPart &gridPart_;
110 template<
class Gr
idPart >
111 template<
class Vector >
112 struct LineSegmentSampler< GridPart >::Reduce
114 Vector
operator() (
const Vector &a,
const Vector &b )
const
117 for(
int k = 0; k < Vector::dimension; ++k )
118 c[ k ] =
std::min( a[ k ], b[ k ] );
128 template<
class Gr
idPart >
129 template<
class Gr
idFunction >
130 inline void LineSegmentSampler< GridPart >
131 ::operator() (
const GridFunction &f, std::vector< typename GridFunction::RangeType > &samples )
const
133 typedef typename GridPartType::template Codim< 0 >::IteratorType IteratorType;
134 typedef typename GridPartType::template Codim< 0 >::EntityType EntityType;
136 const int numSamples = samples.size();
138 DUNE_THROW( InvalidStateException,
"LineSegmentSampler cannot sample less than 2 points." );
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 );
149 const IteratorType end = gridPart().template end< 0 >();
150 for( IteratorType it = gridPart().
template begin< 0 >(); it != end; ++it )
152 const EntityType &entity = *it;
153 const typename EntityType::Geometry &geometry = entity.geometry();
158 auto &refElement = ReferenceElements::general( geometry.type() );
159 const int numFaces = refElement.size( 1 );
160 for(
int face = 0; face < numFaces; ++face )
165 const LocalDomainType &refNormal = refElement.integrationOuterNormal( face );
166 geometry.jacobianInverseTransposed( center ).mv( refNormal, normal );
169 const DomainFieldType ncl = normal * (geometry.global( center ) - left_);
179 else if( ncl < -1e-8 )
182 lower = std::numeric_limits< DomainFieldType >::max();
183 upper = std::numeric_limits< DomainFieldType >::min();
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 )
196 lf.evaluate( geometry.local( xi ), samples[ i ] );
202 typedef Reduce< typename GridFunction::RangeType > Op;
203 gridPart().comm().template allreduce< Op >( &(samples[ 0 ]), numSamples );
208 if constexpr ( std::is_floating_point< RangeFieldType >::value )
210 for(
int i = 0; i < numSamples; ++i )
212 for(
int d=0; d<GridFunction::RangeType::dimension; ++d )
214 valid &= ! std::isnan( samples[ i ][ d ] );
220 DUNE_THROW( InvalidStateException,
"LineSegmentSampler could not find all samples." );
224 template<
class Gr
idPart >
228 const int numSamples = points.size();
230 DUNE_THROW( InvalidStateException,
"LineSegmentSampler cannot sample less than 2 points." );
234 for(
int i = 0; i < numSamples; ++i )
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