Dune Core Modules (2.8.0)

raviartthomassimplexprebasis.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_RAVIARTTHOMASPREBASIS_HH
4#define DUNE_RAVIARTTHOMASPREBASIS_HH
5
6#include <fstream>
7#include <utility>
8
10
11#include <dune/localfunctions/utility/polynomialbasis.hh>
12
13namespace Dune
14{
15 template < GeometryType::Id geometryId, class Field >
16 struct RTVecMatrix;
17
18 template <unsigned int dim, class Field>
19 struct RTPreBasisFactory
20 {
21 typedef MonomialBasisProvider<dim,Field> MBasisFactory;
22 typedef typename MBasisFactory::Object MBasis;
23 typedef StandardEvaluator<MBasis> EvalMBasis;
24 typedef PolynomialBasisWithMatrix<EvalMBasis,SparseCoeffMatrix<Field,dim> > Basis;
25
26 typedef const Basis Object;
27 typedef std::size_t Key;
28
29 template <unsigned int dd, class FF>
30 struct EvaluationBasisFactory
31 {
32 typedef MonomialBasisProvider<dd,FF> Type;
33 };
34 template< GeometryType::Id geometryId >
35 static Object *create ( const Key &order )
36 {
37 RTVecMatrix<geometryId,Field> vecMatrix(order);
38 MBasis *mbasis = MBasisFactory::template create<geometryId>(order+1);
39 typename std::remove_const<Object>::type *tmBasis = new typename std::remove_const<Object>::type(*mbasis);
40 tmBasis->fill(vecMatrix);
41 return tmBasis;
42 }
43 static void release( Object *object ) { delete object; }
44 };
45
46 template <GeometryType::Id geometryId, class Field>
47 struct RTVecMatrix
48 {
49 static constexpr GeometryType geometry = geometryId;
50 static const unsigned int dim = geometry.dim();
51 typedef MultiIndex<dim,Field> MI;
52 typedef MonomialBasis<geometryId,MI> MIBasis;
53 RTVecMatrix(std::size_t order)
54 {
55 /*
56 * Construction of Raviart-Thomas elements in high dimensions see "Mixed Finite Elements in \R^3" by Nedelec, 1980.
57 *
58 * Let $\P_{n,k}$ be the space of polynomials in $n$ variables with degree $\leq k$.
59 * The space of Raviart-Thomas functions in $n$ dimensions with index $k$ is defined as
60 *
61 * \begin{equation*}
62 * RT_k := (\P_{k-1})^n \oplus \widetilde \P_k x
63 * \end{equation*}
64 * with $x=(x_1,x_2,\dots, x_n)$ in $n$ dimensions and $\widetilde \P_k$ the homogeneous polynomials of degree $k$.
65 *
66 * For $RT_k$ holds
67 * \begin{equation*}
68 * (\P_{k-1})^n \subset RT_k \subset (\P_k)^n.
69 * \end{equation*}
70 *
71 * We construct $(\P_k)^n$ and and only use the monomials contained in $RT_k$.
72 *
73 */
74
75 MIBasis basis(order+1);
76 FieldVector< MI, dim > x;
77 /*
78 * Init MultiIndices
79 * x[0]=(1,0,0) x
80 * x[1]=(0,1,0) y
81 * x[2]=(0,0,1) z
82 */
83 for( unsigned int i = 0; i < dim; ++i )
84 x[ i ].set( i, 1 );
85 std::vector< MI > val( basis.size() );
86
87 // val now contains all monomials in $n$ dimensions with degree $\leq order+1$
88 basis.evaluate( x, val );
89
90 col_ = basis.size();
91
92 // get $\dim (\P_{order-1})$
93 unsigned int notHomogen = 0;
94 if (order>0)
95 notHomogen = basis.sizes()[order-1];
96
97 // get $\dim \widetilde (\P_order)$
98 unsigned int homogen = basis.sizes()[order]-notHomogen;
99
100 /*
101 *
102 * The set $RT_k$ is defined as
103 *
104 * \begin{equation}
105 * RT_k := (\P_k)^dim + \widetilde \P_k x with x\in \R^n.
106 * \end{equation}
107 *
108 * The space $\P_k$ is split in $\P_k = \P_{k-1} + \widetilde \P_k$.
109 *
110 * \begin{align}
111 * RT_k &= (\P_{k-1} + \widetilde \P_k)^dim + \widetilde \P_k x with x\in \R^n
112 * &= (\P_{k-1})^n + (\widetilde \P_k)^n + \widetilde \P_k x with x\in \R^n
113 * \end{align}
114 *
115 * Thus $\dim RT_k = n * \dim \P_{k-1} + (n+1)*\dim (\widetilde \P_k)$
116 */
117
118 // row_ = \dim RT_k *dim
119 row_ = (notHomogen*dim+homogen*(dim+1))*dim;
120 mat_ = new Field*[row_];
121 int row = 0;
122
123 /* Assign the correct values for the nonhomogeneous polymonials $p\in (\P_{oder-1})^dim$
124 * A basis function is represented by $dim$ rows.
125 */
126 for (unsigned int i=0; i<notHomogen+homogen; ++i)
127 {
128 for (unsigned int r=0; r<dim; ++r)
129 {
130 for (unsigned int rr=0; rr<dim; ++rr)
131 {
132 // init row to zero
133 mat_[row] = new Field[col_];
134 for (unsigned int j=0; j<col_; ++j)
135 {
136 mat_[row][j] = 0.;
137 }
138 if (r==rr)
139 mat_[row][i] = 1.;
140 ++row;
141 }
142 }
143 }
144
145 /* Assign the correct values for the homogeneous polymonials $p\in RT_k \backslash (\P_{oder-1})^dim$
146 * A basis function is represented by $dim$ rows.
147 */
148 for (unsigned int i=0; i<homogen; ++i)
149 {
150 for (unsigned int r=0; r<dim; ++r)
151 {
152 // init rows to zero
153 mat_[row] = new Field[col_];
154 for (unsigned int j=0; j<col_; ++j)
155 {
156 mat_[row][j] = 0.;
157 }
158
159 unsigned int w;
160 // get a monomial $ p \in \widetilde \P_{order}$
161 MI xval = val[notHomogen+i];
162 xval *= x[r];
163 for (w=homogen+notHomogen; w<val.size(); ++w)
164 {
165 if (val[w] == xval)
166 {
167 mat_[row][w] = 1.;
168 break;
169 }
170 }
171 assert(w<val.size());
172 ++row;
173 }
174 }
175 }
176
177 ~RTVecMatrix()
178 {
179 for (unsigned int i=0; i<rows(); ++i) {
180 delete [] mat_[i];
181 }
182 delete [] mat_;
183 }
184
185 unsigned int cols() const {
186 return col_;
187 }
188
189 unsigned int rows() const {
190 return row_;
191 }
192
193 template <class Vector>
194 void row( const unsigned int row, Vector &vec ) const
195 {
196 const unsigned int N = cols();
197 assert( vec.size() == N );
198 for (unsigned int i=0; i<N; ++i)
199 field_cast(mat_[row][i],vec[i]);
200 }
201 unsigned int row_,col_;
202 Field **mat_;
203 };
204
205
206}
207#endif // DUNE_RAVIARTTHOMASPREBASIS_HH
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:130
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
A unique label for each type of element that can occur in a grid.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)