1 #ifndef DUNE_FEM_OPERATOR_LINEAR_TEST_CHECKLINEAROPERATOR_HH
2 #define DUNE_FEM_OPERATOR_LINEAR_TEST_CHECKLINEAROPERATOR_HH
3 
4 #include <cstddef>
5 #include <cmath>
6 
7 #include <algorithm>
8 #include <utility>
9 #include <vector>
10 
11 #include <dune/common/classname.hh>
12 #include <dune/common/dynvector.hh>
13 #include <dune/common/typetraits.hh>
14 #include <dune/common/typeutilities.hh>
15 
16 #include <dune/fem/common/hybrid.hh>
17 #include <dune/fem/operator/common/stencil.hh>
18 #include <dune/fem/operator/common/temporarylocalmatrix.hh>
19 
20 namespace Dune
21 {
22 
23   namespace Fem
24   {
25 
26     // DiagonalRange
27     // -------------
28 
29     template< class DomainSpace, class RangeSpace >
30     struct DiagonalRange
31     {
32       struct Iterator
33       {
34         typedef typename DomainSpace::IteratorType Iterator1;
35         typedef typename DomainSpace::EntityType Entity1;
36         typedef typename RangeSpace::IteratorType Iterator2;
37         typedef typename RangeSpace::EntityType Entity2;
38 
39         Iterator () = default;
40 
IteratorDune::Fem::DiagonalRange::Iterator41         Iterator ( Iterator1 it1, Iterator2 it2 ) : it_( it1, it2 ) {}
42 
operator ++Dune::Fem::DiagonalRange::Iterator43         Iterator &operator++ () { ++it_.first; ++it_.second; return *this; }
44 
operator *Dune::Fem::DiagonalRange::Iterator45         std::pair< Entity1, Entity2 > operator* () const { return std::make_pair( *it_.first, *it_.second ); }
46 
operator !=Dune::Fem::DiagonalRange::Iterator47         bool operator!= ( const Iterator &other ) const { return (it_ != other.it_); }
operator ==Dune::Fem::DiagonalRange::Iterator48         bool operator== ( const Iterator &other ) const { return (it_ == other.it_); }
49 
50       private:
51         std::pair< Iterator1, Iterator2 > it_;
52       };
53 
DiagonalRangeDune::Fem::DiagonalRange54       DiagonalRange ( const DomainSpace &dSpace, const RangeSpace &rSpace )
55         : dSpace_( dSpace ), rSpace_( rSpace )
56       {}
57 
beginDune::Fem::DiagonalRange58       Iterator begin () const { return Iterator( dSpace_.begin(), rSpace_.begin() ); }
endDune::Fem::DiagonalRange59       Iterator end () const { return Iterator( dSpace_.end(), rSpace_.end() ); }
60 
61     protected:
62       const DomainSpace &dSpace_;
63       const RangeSpace &rSpace_;
64     };
65 
66 
67 
68     // diagonalRange
69     // -------------
70 
71     template< class Space1, class Space2 >
diagonalRange(const Space1 & space1,const Space2 & space2)72     DiagonalRange< Space1, Space2 > diagonalRange ( const Space1 &space1, const Space2 &space2 )
73     {
74       return DiagonalRange< Space1, Space2 >( space1, space2 );
75     }
76 
77 
78 
79     // RangeStencil
80     // ------------
81 
82     template< class DomainSpace, class RangeSpace >
83     struct RangeStencil
84       : public Dune::Fem::Stencil< DomainSpace, RangeSpace >
85     {
86       typedef RangeStencil< DomainSpace, RangeSpace > ThisType;
87       typedef Dune::Fem::Stencil< DomainSpace, RangeSpace > BaseType;
88 
89       template< class Range >
RangeStencilDune::Fem::RangeStencil90       RangeStencil ( const DomainSpace &dSpace, const RangeSpace &rSpace, const Range &range )
91         : BaseType( dSpace, rSpace )
92       {
93         for( const auto &entry : range )
94           BaseType::fill( entry.first, entry.second );
95       }
96     };
97 
98 
99 
100     namespace CheckLinearOperator
101     {
102 
103       template< class LocalMatrix >
verifyPermutation(const LocalMatrix & localMatrix,const std::vector<std::pair<int,int>> & permutation)104       inline void verifyPermutation ( const LocalMatrix &localMatrix, const std::vector< std::pair< int, int > > &permutation )
105       {
106         for( const auto &p : permutation )
107         {
108           if( localMatrix.get( p.first, p.second ) != 1.0 )
109             DUNE_THROW( Dune::NotImplemented, "LocalMatrix not set correctly" );
110         }
111       }
112 
113       template< class LinearOperator, class DomainEntity, class RangeEntity >
114       inline void_t< typename LinearOperator::LocalMatrixType >
verifyLocalMatrixPermutation(const LinearOperator & linOp,const DomainEntity & domainEntity,const RangeEntity & rangeEntity,const std::vector<std::pair<int,int>> & permutation,PriorityTag<1>)115       verifyLocalMatrixPermutation ( const LinearOperator &linOp, const DomainEntity &domainEntity, const RangeEntity &rangeEntity,
116                                      const std::vector< std::pair< int, int > > &permutation, PriorityTag< 1 > )
117       {
118         Dune::Fem::TemporaryLocalMatrix< typename LinearOperator::DomainSpaceType, typename LinearOperator::RangeSpaceType >
119           localMat( linOp.domainSpace(), linOp.rangeSpace() );
120         localMat.bind( domainEntity, rangeEntity );
121         linOp.getLocalMatrix( domainEntity, rangeEntity, localMat );
122         verifyPermutation( localMat, permutation );
123         localMat.unbind();
124       }
125 
126       template< class LinearOperator, class DomainEntity, class RangeEntity >
verifyLocalMatrixPermutation(const LinearOperator & linOp,const DomainEntity & domainEntity,const RangeEntity & rangeEntity,const std::vector<std::pair<int,int>> & permutation,PriorityTag<0>)127       inline void verifyLocalMatrixPermutation ( const LinearOperator &linOp, const DomainEntity &domainEntity, const RangeEntity &rangeEntity,
128                                                  const std::vector< std::pair< int, int > > &permutation, PriorityTag< 0 > )
129       {
130         static bool warn = true;
131         if( warn )
132           std::cout << "Note: " << className( linOp ) << " does not support old-style local matrix interface." << std::endl;
133         warn = false;
134       }
135 
136       template< class LinearOperator, class DomainEntity, class RangeEntity >
verifyLocalMatrixPermutation(const LinearOperator & linOp,const DomainEntity & domainEntity,const RangeEntity & rangeEntity,const std::vector<std::pair<int,int>> & permutation)137       inline void verifyLocalMatrixPermutation ( const LinearOperator &linOp, const DomainEntity &domainEntity, const RangeEntity &rangeEntity,
138                                                 const std::vector< std::pair< int, int > > &permutation )
139       {
140         verifyLocalMatrixPermutation( linOp, domainEntity, rangeEntity, permutation, PriorityTag< 42 >() );
141       }
142 
143     } // namespace CheckLinearOperator
144 
145 
146 
147     // checkLinearOperator
148     // -------------------
149 
150     template< class LinearOperator, class Range >
checkLinearOperator(LinearOperator & linOp,const Range & range,const std::vector<std::pair<int,int>> & permutation)151     inline void checkLinearOperator ( LinearOperator &linOp, const Range &range, const std::vector< std::pair< int, int > > &permutation )
152     {
153       typedef typename LinearOperator::DomainFunctionType DomainFunctionType;
154       typedef typename LinearOperator::RangeFunctionType RangeFunctionType;
155 
156       // get type of domain and range function
157       typedef typename DomainFunctionType::DiscreteFunctionSpaceType DomainSpaceType;
158       typedef typename RangeFunctionType::DiscreteFunctionSpaceType RangeSpaceType;
159 
160       // get spaces
161       const DomainSpaceType &dSpace = linOp.domainSpace();
162       const RangeSpaceType &rSpace = linOp.rangeSpace();
163 
164       // construct some domain and range functions
165       DomainFunctionType domainFunction( "domain Function", dSpace );
166       RangeFunctionType rangeFunction1( "range Function1", rSpace );
167       RangeFunctionType rangeFunction2( "range Function2", rSpace );
168 
169       // first test initialization
170       RangeStencil< DomainSpaceType, RangeSpaceType > stencil( dSpace, rSpace, range );
171 
172       for( int k=0; k<4; ++k )
173       {
174         // check setting operator entries to zero
175         linOp.reserve( stencil );
176         linOp.clear();
177 
178         Dune::Fem::TemporaryLocalMatrix< DomainSpaceType, RangeSpaceType > temp( dSpace, rSpace );
179 
180         // test assemble
181         // check {add,addScaled}LocalMatrix
182         for( const auto &entry : range )
183         {
184           auto domainEntity = entry.first;
185           auto rangeEntity = entry.second;
186 
187           {
188             temp.init( domainEntity, rangeEntity );
189             temp.clear();
190             for( const auto &p : permutation )
191               temp.set( p.first, p.second, 1.0 );
192 
193             linOp.addLocalMatrix( domainEntity, rangeEntity, temp );
194             linOp.addScaledLocalMatrix( domainEntity, rangeEntity, temp, -1.0 );
195           }
196         }
197 
198         linOp.flushAssembly();
199 
200         // check {set}LocalMatrix
201         for( const auto &entry : range )
202         {
203           auto domainEntity = entry.first;
204           auto rangeEntity = entry.second;
205 
206           {
207             // check {add,set,addScaled}LocalMatrix
208             temp.init( domainEntity, rangeEntity );
209             temp.clear();
210             for( const auto &p : permutation )
211               temp.set( p.first, p.second, 1.0 );
212 
213             linOp.setLocalMatrix( domainEntity, rangeEntity, temp );
214           }
215         }
216 
217         linOp.flushAssembly();
218         if (k>=2) // just for checking 'clearRow' and reassembly so don't check content of matrix
219         {
220           try
221           {
222             linOp.unitRow(0);
223             linOp.finalize();
224           }
225           catch( const Dune::NotImplemented &exception )
226           {
227             std::cout << "operator does not implement unitRow method\n";
228           }
229           continue;
230         }
231 
232         // test getLocalMatrix
233         for( const auto &entry : range )
234         {
235           auto domainEntity = entry.first;
236           auto rangeEntity = entry.second;
237 
238           temp.init( domainEntity, rangeEntity );
239           temp.clear();
240           linOp.getLocalMatrix( domainEntity, rangeEntity, temp );
241           CheckLinearOperator::verifyPermutation( temp, permutation );
242 
243           // check old localMatrix implementation
244           CheckLinearOperator::verifyLocalMatrixPermutation( linOp, domainEntity, rangeEntity, permutation );
245         }
246 
247         // test application
248         domainFunction.clear();
249         std::size_t minNumDofs = std::min( domainFunction.blocks(), rangeFunction1.blocks() );
250         for( std::size_t i = 0; i < minNumDofs; ++i )
251         {
252           auto block = domainFunction.dofVector()[ i ];
253           Hybrid::forEach( typename DomainSpaceType::LocalBlockIndices(), [ &block, i ] ( auto &&j ) { block[ j ] = i; } );
254         }
255         linOp( domainFunction, rangeFunction1 );
256 
257         // mimic operator apply manually
258         rangeFunction2.clear();
259         Dune::DynamicVector< double > tmp;
260         std::vector< double > values;
261         for( const auto &entry : range )
262         {
263           auto domainEntity = entry.first;
264           auto rangeEntity = entry.second;
265 
266           tmp.resize( dSpace.blockMapper().numDofs( domainEntity ) * DomainSpaceType::localBlockSize );
267           domainFunction.getLocalDofs( domainEntity, tmp );
268 
269           values.clear();
270           values.resize( rSpace.blockMapper().numDofs( rangeEntity ) * RangeSpaceType::localBlockSize, 0.0 );
271 
272           for( const auto &p : permutation )
273             values[ p.first ] = tmp[ p.second ];
274           rangeFunction2.setLocalDofs( rangeEntity, values );
275         }
276         rangeFunction2.communicate();
277 
278         rangeFunction1.dofVector() -= rangeFunction2.dofVector();
279 
280         double diff = rangeFunction1.scalarProductDofs( rangeFunction1 );
281         if( std::abs( diff ) > 1e-8 )
282           DUNE_THROW( Dune::NotImplemented, "Id Operator not assembled correctly" );
283       }
284     }
285 
286   } // namespace Fem
287 
288 } // namespace Dune
289 
290 #endif // #ifndef DUNE_FEM_OPERATOR_LINEAR_TEST_CHECKLINEAROPERATOR_HH
291