3 #ifndef DUNE_DUAL_Q1_LOCALFINITEELEMENT_HH
4 #define DUNE_DUAL_Q1_LOCALFINITEELEMENT_HH
6 #include <dune/common/fvector.hh>
7 #include <dune/common/fmatrix.hh>
9 #include <dune/geometry/type.hh>
10 #include <dune/geometry/quadraturerules.hh>
25 template<
class D,
class R,
int dim>
43 const Dune::QuadratureRule<D,dim>& quad = Dune::QuadratureRules<D,dim>::rule(gt, 2*dim);
45 const int size = 1 <<dim;
48 Dune::FieldMatrix<R, size, size> massMat;
52 std::vector<Dune::FieldVector<R,1> > integral(size);
53 for (
int i=0; i<size; i++)
56 for(
size_t pt=0; pt<quad.size(); pt++) {
58 const Dune::FieldVector<D ,dim>& pos = quad[pt].position();
59 std::vector<typename Traits::LocalBasisType::Traits::RangeType> q1Values(size);
62 for (
int i=0; i<size; i++) {
66 for (
int j=0; j<dim; j++)
68 q1Values[i] *= (i & (1<<j)) ? pos[j] : 1-pos[j];
71 double weight = quad[pt].weight();
73 for (
int k=0; k<size; k++) {
74 integral[k] += q1Values[k]*weight;
76 for (
int l=0; l<=k; l++)
77 massMat[k][l] += weight*(q1Values[k]*q1Values[l]);
82 for (
int i=0; i<size-1; i++)
83 for (
int j=i+1; j<size; j++)
84 massMat[i][j] = massMat[j][i];
87 Dune::array<Dune::FieldVector<R, size>, size> coefficients;
89 for (
int i=0; i<size; i++) {
91 Dune::FieldVector<R, size> rhs(0);
95 massMat.solve(coefficients[i] ,rhs);
99 basis.setCoefficients(coefficients);
121 return interpolation;