dune-localfunctions  2.4.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 <array>
7 
8 #include <dune/common/fvector.hh>
9 #include <dune/common/fmatrix.hh>
10 
11 #include <dune/geometry/type.hh>
12 #include <dune/geometry/referenceelements.hh>
13 #include <dune/geometry/quadraturerules.hh>
14 
20 
21 namespace Dune
22 {
38  template<class D, class R, int dim, bool faceDual=false>
40  {
41  public:
46 
50  {
51  gt.makeCube(dim);
52 
53  if (faceDual)
54  setupFaceDualCoefficients();
55  else
56  setupDualCoefficients();
57  }
58 
61  const typename Traits::LocalBasisType& localBasis () const
62  {
63  return basis;
64  }
65 
69  {
70  return coefficients;
71  }
72 
76  {
77  return interpolation;
78  }
79 
81  unsigned int size () const
82  {
83  return basis.size();
84  }
85 
88  GeometryType type () const
89  {
90  return gt;
91  }
92 
94  {
95  return new DualQ1LocalFiniteElement(*this);
96  }
97 
98  private:
100  void setupFaceDualCoefficients();
101 
103  void setupDualCoefficients();
104 
106  DualQ1LocalCoefficients<dim> coefficients;
108  GeometryType gt;
109  };
110 
111  template<class D, class R, int dim, bool faceDual>
113  {
114 
115  const int size = 1 <<dim;
116  std::array<Dune::FieldVector<R, size>, size> coeffs;
117 
118  // dual basis functions are linear combinations of Lagrange elements
119  // compute these coefficients here because the basis and the local interpolation needs them
120  const Dune::QuadratureRule<D,dim>& quad = Dune::QuadratureRules<D,dim>::rule(gt, 2*dim);
121 
122  // assemble mass matrix on the reference element
123  Dune::FieldMatrix<R, size, size> massMat;
124  massMat = 0;
125 
126  // and the integrals of the lagrange shape functions
127  std::vector<Dune::FieldVector<R,1> > integral(size);
128  for (int i=0; i<size; i++)
129  integral[i] = 0;
130 
132  for(size_t pt=0; pt<quad.size(); pt++) {
133 
134  const Dune::FieldVector<D ,dim>& pos = quad[pt].position();
135  std::vector<Dune::FieldVector<R,1> > q1Values(size);
136  q1Basis.evaluateFunction(pos,q1Values);
137 
138  D weight = quad[pt].weight();
139 
140  for (int k=0; k<size; k++) {
141  integral[k] += q1Values[k]*weight;
142 
143  for (int l=0; l<=k; l++)
144  massMat[k][l] += weight*(q1Values[k]*q1Values[l]);
145  }
146  }
147 
148  // make matrix symmetric
149  for (int i=0; i<size-1; i++)
150  for (int j=i+1; j<size; j++)
151  massMat[i][j] = massMat[j][i];
152 
153  //solve for the coefficients
154 
155  for (int i=0; i<size; i++) {
156 
157  Dune::FieldVector<R, size> rhs(0);
158  rhs[i] = integral[i];
159 
160  coeffs[i] = 0;
161  massMat.solve(coeffs[i] ,rhs);
162 
163  }
164 
165  basis.setCoefficients(coeffs);
166  interpolation.setCoefficients(coeffs);
167  }
168 
169  template<class D, class R, int dim, bool faceDual>
171  {
172 
173  const int size = 1 <<dim;
174  std::array<Dune::FieldVector<R, size>, size> coeffs;
175 
176  // dual basis functions are linear combinations of Lagrange elements
178 
179  typedef Dune::ReferenceElement<D,dim> RefElement;
180  const RefElement& refElement = Dune::ReferenceElements<D,dim>::general(gt);
181 
182  // loop over faces
183  for (size_t i=0; i<refElement.size(1);i++) {
184 
185  const Dune::QuadratureRule<D,dim-1>& quad = Dune::QuadratureRules<D,dim-1>::rule(refElement.type(i,1),2*dim);
186 
187  // for each face assemble the mass matrix over the face of all
188  // non-vanishing basis functions,
189  // for cubes refElement.size(i,1,dim)=size/2
190  Dune::FieldMatrix<R, size/2, size/2> massMat;
191  massMat = 0;
192 
193  // get geometry
194  const typename RefElement::template Codim<1>::Geometry& geometry = refElement.template geometry<1>(i);
195 
196  // and the integrals of the lagrange shape functions
197  std::vector<Dune::FieldVector<R,1> > integral(size/2);
198  for (int k=0; k<size/2; k++)
199  integral[k] = 0;
200 
201  for(size_t pt=0; pt<quad.size(); pt++) {
202 
203  const Dune::FieldVector<D,dim-1>& pos = quad[pt].position();
204  const Dune::FieldVector<D,dim> elementPos = geometry.global(pos);
205 
206  std::vector<Dune::FieldVector<R,1> > q1Values(size);
207  q1Basis.evaluateFunction(elementPos,q1Values);
208 
209  D weight = quad[pt].weight();
210 
211  for (int k=0; k<refElement.size(i,1,dim); k++) {
212  int row = refElement.subEntity(i,1,k,dim);
213  integral[k] += q1Values[row]*weight;
214 
215  for (int l=0; l<refElement.size(i,1,dim); l++) {
216  int col = refElement.subEntity(i,1,l,dim);
217  massMat[k][l] += weight*(q1Values[row]*q1Values[col]);
218  }
219  }
220  }
221 
222  // solve for the coefficients
223  // note that we possibly overwrite coefficients for neighbouring faces
224  // this is okay since the coefficients are symmetric
225  for (int l=0; l<refElement.size(i,1,dim); l++) {
226 
227  int row = refElement.subEntity(i,1,l,dim);
228  Dune::FieldVector<R, size/2> rhs(0);
229  rhs[l] = integral[l];
230 
231  Dune::FieldVector<R, size/2> x(0);
232  massMat.solve(x ,rhs);
233 
234  for (int k=0; k<refElement.size(i,1,dim); k++) {
235  int col = refElement.subEntity(i,1,k,dim);
236  coeffs[row][col]=x[k];
237  }
238  }
239  }
240 
241  basis.setCoefficients(coeffs);
242  interpolation.setCoefficients(coeffs);
243  }
244 }
245 #endif
traits helper struct
Definition: localfiniteelementtraits.hh:10
Layout map for dual Q1 elements.
Definition: dualq1localcoefficients.hh:22
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: dualq1.hh:68
unsigned int size() const
Number of shape functions in this finite element.
Definition: dualq1.hh:81
DualQ1LocalFiniteElement * clone() const
Definition: dualq1.hh:93
Dual Lagrange shape functions of order 1 on the reference cube.
Definition: dualq1localbasis.hh:25
const Traits::LocalBasisType & localBasis() const
Definition: dualq1.hh:61
const Traits::LocalInterpolationType & localInterpolation() const
Definition: dualq1.hh:75
Definition: brezzidouglasmarini1cube2dlocalbasis.hh:14
LI LocalInterpolationType
Definition: localfiniteelementtraits.hh:22
LC LocalCoefficientsType
Definition: localfiniteelementtraits.hh:18
Definition: dualq1localinterpolation.hh:17
GeometryType type() const
Definition: dualq1.hh:88
Lagrange shape functions of order 1 on the reference cube.
Definition: q1localbasis.hh:24
LB LocalBasisType
Definition: localfiniteelementtraits.hh:14
The local dual Q1 finite element on cubes.
Definition: dualq1.hh:39
LocalFiniteElementTraits< DualQ1LocalBasis< D, R, dim >, DualQ1LocalCoefficients< dim >, DualQ1LocalInterpolation< dim, DualQ1LocalBasis< D, R, dim > > > Traits
Definition: dualq1.hh:45
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: q1localbasis.hh:37
DualQ1LocalFiniteElement()
Definition: dualq1.hh:49