5#ifndef DUNE_LOCALFUNCTIONS_NEDELEC_NEDELECSIMPLEX_NEDELECSIMPLEXPREBASIS_HH
6#define DUNE_LOCALFUNCTIONS_NEDELEC_NEDELECSIMPLEX_NEDELECSIMPLEXPREBASIS_HH
13#include <dune/localfunctions/utility/polynomialbasis.hh>
17 template < GeometryType::Id geometryId,
class Field >
18 struct NedelecVecMatrix;
23 template <
unsigned int dim,
class Field>
24 struct NedelecPreBasisFactory
26 typedef MonomialBasisProvider<dim,Field> MBasisFactory;
27 typedef typename MBasisFactory::Object MBasis;
28 typedef StandardEvaluator<MBasis> EvalMBasis;
29 typedef PolynomialBasisWithMatrix<EvalMBasis,SparseCoeffMatrix<Field,dim> > Basis;
31 typedef const Basis Object;
32 typedef std::size_t Key;
34 template <
unsigned int dd,
class FF>
35 struct EvaluationBasisFactory
37 typedef MonomialBasisProvider<dd,FF> Type;
40 template< GeometryType::Id geometryId >
41 static Object *create ( Key order )
53 NedelecVecMatrix<geometryId,Field> vecMatrix(order);
54 MBasis *mbasis = MBasisFactory::template create<geometryId>(order+1);
55 std::remove_const_t<Object>* tmBasis =
new std::remove_const_t<Object>(*mbasis);
56 tmBasis->fill(vecMatrix);
59 static void release( Object *
object ) {
delete object; }
62 template <GeometryType::Id geometryId,
class Field>
63 struct NedelecVecMatrix
66 static const unsigned int dim = geometry.dim();
67 typedef MultiIndex<dim,Field> MI;
68 typedef MonomialBasis<geometryId,MI> MIBasis;
69 NedelecVecMatrix(std::size_t order)
90 if( (dim!=2 && dim!=3) || !geometry.isSimplex())
93 MIBasis basis(order+1);
94 FieldVector< MI, dim > x;
101 for(
unsigned int i = 0; i < dim; ++i )
103 std::vector< MI > val( basis.size() );
106 basis.evaluate( x, val );
111 unsigned int notHomogen = 0;
113 notHomogen = basis.sizes()[order-1];
116 unsigned int homogen = basis.sizes()[order]-notHomogen;
144 row_ = (notHomogen*dim+homogen*(dim+1))*dim;
148 int homogenTwoVariables = 0;
149 for(
int w = notHomogen; w<notHomogen + homogen; w++)
151 homogenTwoVariables++;
152 row_ = (notHomogen*dim+homogen*(dim+2) + homogenTwoVariables)*dim;
155 mat_ =
new Field*[row_];
161 for (
unsigned int i=0; i<notHomogen+homogen; ++i)
163 for (
unsigned int r=0; r<dim; ++r)
165 for (
unsigned int rr=0; rr<dim; ++rr)
168 mat_[row] =
new Field[col_];
169 for (
unsigned int j=0; j<col_; ++j)
182 for (
unsigned int i=0; i<homogen; ++i)
185 MI xval = val[notHomogen+i];
188 for (
unsigned int r=0; r<dim; ++r)
191 mat_[row+r] =
new Field[col_];
192 for (
unsigned int j=0; j<col_; ++j)
200 for (
int w=homogen+notHomogen; w<val.size(); ++w)
202 if (val[w] == xval*x[0])
204 if (val[w] == xval*x[1])
211 int skipLastDim = xval.z(0)>0;
217 for (
unsigned int r=0; r<dim*(dim-skipLastDim); ++r)
220 mat_[row+r] =
new Field[col_];
221 for (
unsigned int j=0; j<col_; ++j)
233 for (
unsigned int r=0; r<dim - skipLastDim; ++r)
235 int index = (r+dim-1)%dim;
236 for (
int w=homogen+notHomogen; w<val.size(); ++w)
238 if (val[w] == xval*x[index])
240 if (val[w] == xval*x[r])
241 mat_[row+index][w] = -1.;
252 for (
unsigned int i=0; i<rows(); ++i) {
258 unsigned int cols()
const {
262 unsigned int rows()
const {
266 template <
class Vector>
267 void row(
const unsigned int row, Vector &vec )
const
269 const unsigned int N = cols();
270 assert( vec.size() == N );
271 for (
unsigned int i=0; i<N; ++i)
275 unsigned int row_,col_;
Default exception for dummy implementations.
Definition: exceptions.hh:263
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:132
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
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.