dune-localfunctions  2.3.1
dualq1.hh
Go to the documentation of this file.
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_DUAL_Q1_LOCALFINITEELEMENT_HH
4 #define DUNE_DUAL_Q1_LOCALFINITEELEMENT_HH
5 
6 #include <dune/common/fvector.hh>
7 #include <dune/common/fmatrix.hh>
8 
9 #include <dune/geometry/type.hh>
10 #include <dune/geometry/quadraturerules.hh>
11 
16 
17 namespace Dune
18 {
19 
25  template<class D, class R, int dim>
27  {
28  public:
33 
37  {
38  gt.makeCube(dim);
39 
40  // dual basis functions are linear combinations of Lagrange elements
41  // compute these coefficients here because the basis and the local interpolation needs them
42 
43  const Dune::QuadratureRule<D,dim>& quad = Dune::QuadratureRules<D,dim>::rule(gt, 2*dim);
44 
45  const int size = 1 <<dim;
46 
47  // assemble mass matrix on the reference element
48  Dune::FieldMatrix<R, size, size> massMat;
49  massMat = 0;
50 
51  // and the integrals of the lagrange shape functions
52  std::vector<Dune::FieldVector<R,1> > integral(size);
53  for (int i=0; i<size; i++)
54  integral[i] = 0;
55 
56  for(size_t pt=0; pt<quad.size(); pt++) {
57 
58  const Dune::FieldVector<D ,dim>& pos = quad[pt].position();
59  std::vector<typename Traits::LocalBasisType::Traits::RangeType> q1Values(size);
60 
61  // evaluate q1 basis functions
62  for (int i=0; i<size; i++) {
63 
64  q1Values[i] = 1;
65 
66  for (int j=0; j<dim; j++)
67  // if j-th bit of i is set multiply with in[j], else with 1-in[j]
68  q1Values[i] *= (i & (1<<j)) ? pos[j] : 1-pos[j];
69  }
70 
71  double weight = quad[pt].weight();
72 
73  for (int k=0; k<size; k++) {
74  integral[k] += q1Values[k]*weight;
75 
76  for (int l=0; l<=k; l++)
77  massMat[k][l] += weight*(q1Values[k]*q1Values[l]);
78  }
79  }
80 
81  // make matrix symmetric
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];
85 
86  //solve for the coefficients
87  Dune::array<Dune::FieldVector<R, size>, size> coefficients;
88 
89  for (int i=0; i<size; i++) {
90 
91  Dune::FieldVector<R, size> rhs(0);
92  rhs[i] = integral[i];
93 
94  coefficients[i] = 0;
95  massMat.solve(coefficients[i] ,rhs);
96 
97  }
98 
99  basis.setCoefficients(coefficients);
100  interpolation.setCoefficients(coefficients);
101  }
102 
105  const typename Traits::LocalBasisType& localBasis () const
106  {
107  return basis;
108  }
109 
113  {
114  return coefficients;
115  }
116 
120  {
121  return interpolation;
122  }
123 
126  GeometryType type () const
127  {
128  return gt;
129  }
130 
132  {
133  return new DualQ1LocalFiniteElement(*this);
134  }
135 
136  private:
138  DualQ1LocalCoefficients<dim> coefficients;
140  GeometryType gt;
141  };
142 }
143 
144 #endif
traits helper struct
Definition: localfiniteelementtraits.hh:10
DualQ1LocalFiniteElement()
Definition: dualq1.hh:36
Layout map for dual Q1 elements.
Definition: dualq1localcoefficients.hh:22
const Traits::LocalInterpolationType & localInterpolation() const
Definition: dualq1.hh:119
DualQ1LocalFiniteElement * clone() const
Definition: dualq1.hh:131
Dual Lagrange shape functions of order 1 on the reference cube.
Definition: dualq1localbasis.hh:23
const Traits::LocalBasisType & localBasis() const
Definition: dualq1.hh:105
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: dualq1.hh:112
LI LocalInterpolationType
Definition: localfiniteelementtraits.hh:22
LocalFiniteElementTraits< DualQ1LocalBasis< D, R, dim >, DualQ1LocalCoefficients< dim >, DualQ1LocalInterpolation< dim, DualQ1LocalBasis< D, R, dim > > > Traits
Definition: dualq1.hh:32
LC LocalCoefficientsType
Definition: localfiniteelementtraits.hh:18
Definition: dualq1localinterpolation.hh:16
GeometryType type() const
Definition: dualq1.hh:126
LB LocalBasisType
Definition: localfiniteelementtraits.hh:14
The local dual Q1 finite element on cubes.
Definition: dualq1.hh:26