DUNE-FEM (unstable)

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