DUNE PDELab (git)

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// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
4// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5#ifndef DUNE_BASISMATRIX_HH
6#define DUNE_BASISMATRIX_HH
7
8#include <fstream>
11
12#include <dune/localfunctions/utility/monomialbasis.hh>
13#include <dune/localfunctions/utility/polynomialbasis.hh>
14
15namespace Dune
16{
17 /****************************************
18 * A dense matrix representation of a ''polynomial''
19 * basis. Its represent a basis as a linear
20 * combination of a second basis, i.e., a
21 * monomial basis. It is simular to the PolynomialBasis
22 * but it not derived from the LocalBasis class.
23 * It is used to define a ''pre basis''.
24 ****************************************/
25 template< class PreBasis, class Interpolation,
26 class Field >
27 struct BasisMatrix;
28
29 template< class PreBasis, class Interpolation,
30 class Field >
31 struct BasisMatrixBase : public DynamicMatrix<Field>
32 {
33 typedef DynamicMatrix<Field> Matrix;
34
35 BasisMatrixBase( const PreBasis& preBasis,
36 const Interpolation& localInterpolation )
37 : cols_(preBasis.size())
38 {
39 localInterpolation.interpolate( preBasis, *this );
40 this->invert();
41 }
42 unsigned int cols () const
43 {
44 return cols_;
45 }
46 unsigned int rows () const
47 {
48 return Matrix::rows();
49 }
50 private:
51 unsigned int cols_;
52 };
53
54 template< GeometryType::Id geometryId, class F,
55 class Interpolation,
56 class Field >
57 struct BasisMatrix< const MonomialBasis< geometryId, F >, Interpolation, Field >
58 : public BasisMatrixBase< const MonomialBasis< geometryId, F >, Interpolation, Field >
59 {
60 typedef const MonomialBasis< geometryId, F > PreBasis;
61 typedef BasisMatrixBase<PreBasis,Interpolation,Field> Base;
62 typedef typename Base::Matrix Matrix;
63
64 BasisMatrix( const PreBasis& preBasis,
65 const Interpolation& localInterpolation )
66 : Base(preBasis, localInterpolation)
67 {}
68 template <class Vector>
69 void row( const unsigned int row, Vector &vec ) const
70 {
71 const unsigned int N = Matrix::rows();
72 assert( Matrix::cols() == N && vec.size() == N );
73 // note: that the transposed matrix is computed,
74 // and is square
75 for (unsigned int i=0; i<N; ++i)
76 field_cast((*this)[i][row],vec[i]);
77 }
78 };
79 template< int dim, class F,
80 class Interpolation,
81 class Field >
82 struct BasisMatrix< const Dune::VirtualMonomialBasis< dim, F >, Interpolation, Field >
83 : public BasisMatrixBase< const VirtualMonomialBasis< dim, F >, Interpolation, Field >
84 {
85 typedef const VirtualMonomialBasis< dim, F > PreBasis;
86 typedef BasisMatrixBase<PreBasis,Interpolation,Field> Base;
87 typedef typename Base::Matrix Matrix;
88
89 BasisMatrix( const PreBasis& preBasis,
90 const Interpolation& localInterpolation )
91 : Base(preBasis, localInterpolation)
92 {}
93 template <class Vector>
94 void row( const unsigned int row, Vector &vec ) const
95 {
96 const unsigned int N = Matrix::rows();
97 assert( Matrix::cols() == N && vec.size() == N );
98 // note: that the transposed matrix is computed,
99 // and is square
100 for (unsigned int i=0; i<N; ++i)
101 field_cast((*this)[i][row],vec[i]);
102 }
103 };
104 template< class Eval, class CM, class D, class R,
105 class Interpolation,
106 class Field >
107 struct BasisMatrix< const PolynomialBasis<Eval,CM,D,R>, Interpolation, Field >
108 : public BasisMatrixBase< const PolynomialBasis<Eval,CM,D,R>, Interpolation, Field >
109 {
110 typedef const PolynomialBasis<Eval,CM,D,R> PreBasis;
111 typedef BasisMatrixBase<PreBasis,Interpolation,Field> Base;
112 typedef typename Base::Matrix Matrix;
113
114 BasisMatrix( const PreBasis& preBasis,
115 const Interpolation& localInterpolation )
116 : Base(preBasis, localInterpolation),
117 preBasis_(preBasis)
118 {}
119 unsigned int cols() const
120 {
121 return preBasis_.matrix().baseSize() ;
122 }
123 template <class Vector>
124 void row( const unsigned int row, Vector &vec ) const
125 {
126 assert( Matrix::rows() == Matrix::cols() );
127 assert( vec.size() == preBasis_.matrix().baseSize() );
128 assert( Matrix::cols() == preBasis_.size() );
129 for (unsigned int j=0; j<Matrix::cols(); ++j)
130 vec[j] = 0;
131 for (unsigned int i=0; i<Matrix::rows(); ++i)
132 preBasis_.matrix().
133 addRow(i,(*this)[i][row],vec);
134 }
135 private:
136 const PreBasis& preBasis_;
137 };
138 template< class Eval, class CM,
139 class Interpolation,
140 class Field >
141 struct BasisMatrix< const PolynomialBasisWithMatrix<Eval,CM>, Interpolation, Field >
142 : public BasisMatrixBase< const PolynomialBasisWithMatrix<Eval,CM>, Interpolation, Field >
143 {
144 typedef const PolynomialBasisWithMatrix<Eval,CM> PreBasis;
145 typedef BasisMatrixBase<PreBasis,Interpolation,Field> Base;
146 typedef typename Base::Matrix Matrix;
147
148 BasisMatrix( const PreBasis& preBasis,
149 const Interpolation& localInterpolation )
150 : Base(preBasis, localInterpolation),
151 preBasis_(preBasis)
152 {}
153 unsigned int cols() const
154 {
155 return preBasis_.matrix().baseSize() ;
156 }
157 unsigned int rows () const
158 {
159 assert( Matrix::rows() == preBasis_.matrix().size() );
160 return preBasis_.matrix().size()*CM::blockSize ;
161 }
162 template <class Vector>
163 void row( const unsigned int row, Vector &vec ) const
164 {
165 unsigned int r = row / CM::blockSize;
166 assert( r < Matrix::rows() );
167 assert( Matrix::rows() == Matrix::cols() );
168 assert( vec.size() == preBasis_.matrix().baseSize() );
169 assert( Matrix::cols() == preBasis_.size() );
170 for (unsigned int j=0; j<vec.size(); ++j)
171 vec[j] = 0;
172 for (unsigned int i=0; i<Matrix::rows(); ++i)
173 preBasis_.matrix().
174 addRow(i*CM::blockSize+row%CM::blockSize,(*this)[i][r],vec);
175 }
176 private:
177 const PreBasis& preBasis_;
178 };
179}
180
181#endif // DUNE_BASISMATRIX_HH
void invert(bool doPivoting=true)
Compute inverse.
constexpr size_type rows() const
number of rows
Definition: densematrix.hh:709
IdType Id
An integral id representing a GeometryType.
Definition: type.hh:181
A few common exception classes.
This file implements a dense matrix with dynamic numbers of rows and columns.
Dune namespace.
Definition: alignedallocator.hh:13
void field_cast(const F1 &f1, F2 &f2)
a helper class to cast from one field to another
Definition: field.hh:159
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)