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