dune-localfunctions 3.0-git
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
21namespace 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>
112 void DualQ1LocalFiniteElement<D,R,dim,faceDual>::setupDualCoefficients()
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 auto& 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>
170 void DualQ1LocalFiniteElement<D,R,dim,faceDual>::setupFaceDualCoefficients()
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 const auto& refElement = Dune::ReferenceElements<D,dim>::general(gt);
180
181 // loop over faces
182 for (size_t i=0; i<refElement.size(1);i++) {
183
184 const auto& quad = Dune::QuadratureRules<D,dim-1>::rule(refElement.type(i,1),2*dim);
185
186 // for each face assemble the mass matrix over the face of all
187 // non-vanishing basis functions,
188 // for cubes refElement.size(i,1,dim)=size/2
189 Dune::FieldMatrix<R, size/2, size/2> massMat;
190 massMat = 0;
191
192 // get geometry
193 const auto& geometry = refElement.template geometry<1>(i);
194
195 // and the integrals of the lagrange shape functions
196 std::vector<Dune::FieldVector<R,1> > integral(size/2);
197 for (int k=0; k<size/2; k++)
198 integral[k] = 0;
199
200 for(size_t pt=0; pt<quad.size(); pt++) {
201
202 const auto& pos = quad[pt].position();
203 const auto& elementPos = geometry.global(pos);
204
205 std::vector<Dune::FieldVector<R,1> > q1Values(size);
206 q1Basis.evaluateFunction(elementPos,q1Values);
207
208 D weight = quad[pt].weight();
209
210 for (int k=0; k<refElement.size(i,1,dim); k++) {
211 int row = refElement.subEntity(i,1,k,dim);
212 integral[k] += q1Values[row]*weight;
213
214 for (int l=0; l<refElement.size(i,1,dim); l++) {
215 int col = refElement.subEntity(i,1,l,dim);
216 massMat[k][l] += weight*(q1Values[row]*q1Values[col]);
217 }
218 }
219 }
220
221 // solve for the coefficients
222 // note that we possibly overwrite coefficients for neighbouring faces
223 // this is okay since the coefficients are symmetric
224 for (int l=0; l<refElement.size(i,1,dim); l++) {
225
226 int row = refElement.subEntity(i,1,l,dim);
227 Dune::FieldVector<R, size/2> rhs(0);
228 rhs[l] = integral[l];
229
230 Dune::FieldVector<R, size/2> x(0);
231 massMat.solve(x ,rhs);
232
233 for (int k=0; k<refElement.size(i,1,dim); k++) {
234 int col = refElement.subEntity(i,1,k,dim);
235 coeffs[row][col]=x[k];
236 }
237 }
238 }
239
240 basis.setCoefficients(coeffs);
241 interpolation.setCoefficients(coeffs);
242 }
243}
244#endif
Definition brezzidouglasmarini1cube2d.hh:14
traits helper struct
Definition localfiniteelementtraits.hh:11
LB LocalBasisType
Definition localfiniteelementtraits.hh:14
LC LocalCoefficientsType
Definition localfiniteelementtraits.hh:18
LI LocalInterpolationType
Definition localfiniteelementtraits.hh:22
The local dual Q1 finite element on cubes.
Definition dualq1.hh:40
DualQ1LocalFiniteElement * clone() const
Definition dualq1.hh:93
LocalFiniteElementTraits< DualQ1LocalBasis< D, R, dim >, DualQ1LocalCoefficients< dim >, DualQ1LocalInterpolation< dim, DualQ1LocalBasis< D, R, dim > > > Traits
Definition dualq1.hh:45
unsigned int size() const
Number of shape functions in this finite element.
Definition dualq1.hh:81
const Traits::LocalInterpolationType & localInterpolation() const
Definition dualq1.hh:75
const Traits::LocalBasisType & localBasis() const
Definition dualq1.hh:61
const Traits::LocalCoefficientsType & localCoefficients() const
Definition dualq1.hh:68
DualQ1LocalFiniteElement()
Definition dualq1.hh:49
GeometryType type() const
Definition dualq1.hh:88
Dual Lagrange shape functions of order 1 on the reference cube.
Definition dualq1localbasis.hh:27
Layout map for dual Q1 elements.
Definition dualq1localcoefficients.hh:23
Definition dualq1localinterpolation.hh:18
Lagrange shape functions of order 1 on the reference cube.
Definition q1localbasis.hh:27
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition q1localbasis.hh:39