Dune Core Modules (2.8.0)

basismatrix.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_BASISMATRIX_HH
4#define DUNE_BASISMATRIX_HH
5
6#include <fstream>
8
9#include <dune/localfunctions/utility/lfematrix.hh>
10#include <dune/localfunctions/utility/monomialbasis.hh>
11#include <dune/localfunctions/utility/polynomialbasis.hh>
12
13namespace Dune
14{
15 /****************************************
16 * A dense matrix representation of a ''polynomial''
17 * basis. Its represent a basis as a linear
18 * combination of a second basis, i.e., a
19 * monomial basis. It is simular to the PolynomialBasis
20 * but it not derived from the LocalBasis class.
21 * It is used to define a ''pre basis''.
22 ****************************************/
23 template< class PreBasis, class Interpolation,
24 class Field >
25 struct BasisMatrix;
26
27 template< class PreBasis, class Interpolation,
28 class Field >
29 struct BasisMatrixBase : public LFEMatrix<Field>
30 {
31 typedef LFEMatrix<Field> Matrix;
32
33 BasisMatrixBase( const PreBasis& preBasis,
34 const Interpolation& localInterpolation )
35 : cols_(preBasis.size())
36 {
37 localInterpolation.interpolate( preBasis, *this );
38
39 if ( !Matrix::invert() )
40 {
41 DUNE_THROW(MathError, "While computing basis a singular matrix was constructed!");
42 }
43 }
44 unsigned int cols () const
45 {
46 return cols_;
47 }
48 unsigned int rows () const
49 {
50 return Matrix::rows();
51 }
52 private:
53 unsigned int cols_;
54 };
55
56 template< GeometryType::Id geometryId, class F,
57 class Interpolation,
58 class Field >
59 struct BasisMatrix< const MonomialBasis< geometryId, F >, Interpolation, Field >
60 : public BasisMatrixBase< const MonomialBasis< geometryId, F >, Interpolation, Field >
61 {
62 typedef const MonomialBasis< geometryId, F > PreBasis;
63 typedef BasisMatrixBase<PreBasis,Interpolation,Field> Base;
64 typedef typename Base::Matrix Matrix;
65
66 BasisMatrix( const PreBasis& preBasis,
67 const Interpolation& localInterpolation )
68 : Base(preBasis, localInterpolation)
69 {}
70 template <class Vector>
71 void row( const unsigned int row, Vector &vec ) const
72 {
73 const unsigned int N = Matrix::rows();
74 assert( Matrix::cols() == N && vec.size() == N );
75 // note: that the transposed matrix is computed,
76 // and is square
77 for (unsigned int i=0; i<N; ++i)
78 field_cast(Matrix::operator()(i,row),vec[i]);
79 }
80 };
81 template< int dim, class F,
82 class Interpolation,
83 class Field >
84 struct BasisMatrix< const Dune::VirtualMonomialBasis< dim, F >, Interpolation, Field >
85 : public BasisMatrixBase< const VirtualMonomialBasis< dim, F >, Interpolation, Field >
86 {
87 typedef const VirtualMonomialBasis< dim, F > PreBasis;
88 typedef BasisMatrixBase<PreBasis,Interpolation,Field> Base;
89 typedef typename Base::Matrix Matrix;
90
91 BasisMatrix( const PreBasis& preBasis,
92 const Interpolation& localInterpolation )
93 : Base(preBasis, localInterpolation)
94 {}
95 template <class Vector>
96 void row( const unsigned int row, Vector &vec ) const
97 {
98 const unsigned int N = Matrix::rows();
99 assert( Matrix::cols() == N && vec.size() == N );
100 // note: that the transposed matrix is computed,
101 // and is square
102 for (unsigned int i=0; i<N; ++i)
103 field_cast(Matrix::operator()(i,row),vec[i]);
104 }
105 };
106 template< class Eval, class CM, class D, class R,
107 class Interpolation,
108 class Field >
109 struct BasisMatrix< const PolynomialBasis<Eval,CM,D,R>, Interpolation, Field >
110 : public BasisMatrixBase< const PolynomialBasis<Eval,CM,D,R>, Interpolation, Field >
111 {
112 typedef const PolynomialBasis<Eval,CM,D,R> PreBasis;
113 typedef BasisMatrixBase<PreBasis,Interpolation,Field> Base;
114 typedef typename Base::Matrix Matrix;
115
116 BasisMatrix( const PreBasis& preBasis,
117 const Interpolation& localInterpolation )
118 : Base(preBasis, localInterpolation),
119 preBasis_(preBasis)
120 {}
121 unsigned int cols() const
122 {
123 return preBasis_.matrix().baseSize() ;
124 }
125 template <class Vector>
126 void row( const unsigned int row, Vector &vec ) const
127 {
128 assert( Matrix::rows() == Matrix::cols() );
129 assert( vec.size() == preBasis_.matrix().baseSize() );
130 assert( Matrix::cols() == preBasis_.size() );
131 for (unsigned int j=0; j<Matrix::cols(); ++j)
132 vec[j] = 0;
133 for (unsigned int i=0; i<Matrix::rows(); ++i)
134 preBasis_.matrix().
135 addRow(i,Base::Matrix::operator()(i,row),vec);
136 }
137 private:
138 const PreBasis& preBasis_;
139 };
140 template< class Eval, class CM,
141 class Interpolation,
142 class Field >
143 struct BasisMatrix< const PolynomialBasisWithMatrix<Eval,CM>, Interpolation, Field >
144 : public BasisMatrixBase< const PolynomialBasisWithMatrix<Eval,CM>, Interpolation, Field >
145 {
146 typedef const PolynomialBasisWithMatrix<Eval,CM> PreBasis;
147 typedef BasisMatrixBase<PreBasis,Interpolation,Field> Base;
148 typedef typename Base::Matrix Matrix;
149
150 BasisMatrix( const PreBasis& preBasis,
151 const Interpolation& localInterpolation )
152 : Base(preBasis, localInterpolation),
153 preBasis_(preBasis)
154 {}
155 unsigned int cols() const
156 {
157 return preBasis_.matrix().baseSize() ;
158 }
159 unsigned int rows () const
160 {
161 assert( Matrix::rows() == preBasis_.matrix().size() );
162 return preBasis_.matrix().size()*CM::blockSize ;
163 }
164 template <class Vector>
165 void row( const unsigned int row, Vector &vec ) const
166 {
167 unsigned int r = row / CM::blockSize;
168 assert( r < Matrix::rows() );
169 assert( Matrix::rows() == Matrix::cols() );
170 assert( vec.size() == preBasis_.matrix().baseSize() );
171 assert( Matrix::cols() == preBasis_.size() );
172 for (unsigned int j=0; j<vec.size(); ++j)
173 vec[j] = 0;
174 for (unsigned int i=0; i<Matrix::rows(); ++i)
175 preBasis_.matrix().
176 addRow(i*CM::blockSize+row%CM::blockSize,Base::Matrix::operator()(i,r),vec);
177 }
178 private:
179 const PreBasis& preBasis_;
180 };
181}
182
183#endif // DUNE_BASISMATRIX_HH
IdType Id
An integral id representing a GeometryType.
Definition: type.hh:190
A few common exception classes.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Dune namespace.
Definition: alignedallocator.hh:11
void field_cast(const F1 &f1, F2 &f2)
a helper class to cast from one field to another
Definition: field.hh:157
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)