3 #ifndef DUNE_DUAL_Q1_LOCALFINITEELEMENT_HH
4 #define DUNE_DUAL_Q1_LOCALFINITEELEMENT_HH
8 #include <dune/common/fvector.hh>
9 #include <dune/common/fmatrix.hh>
11 #include <dune/geometry/type.hh>
12 #include <dune/geometry/referenceelements.hh>
13 #include <dune/geometry/quadraturerules.hh>
38 template<
class D,
class R,
int dim,
bool faceDual=false>
54 setupFaceDualCoefficients();
56 setupDualCoefficients();
100 void setupFaceDualCoefficients();
103 void setupDualCoefficients();
106 DualQ1LocalCoefficients<dim> coefficients;
111 template<
class D,
class R,
int dim,
bool faceDual>
112 void DualQ1LocalFiniteElement<D,R,dim,faceDual>::setupDualCoefficients()
115 const int size = 1 <<dim;
116 std::array<Dune::FieldVector<R, size>, size> coeffs;
120 const Dune::QuadratureRule<D,dim>& quad = Dune::QuadratureRules<D,dim>::rule(gt, 2*dim);
123 Dune::FieldMatrix<R, size, size> massMat;
127 std::vector<Dune::FieldVector<R,1> > integral(size);
128 for (
int i=0; i<size; i++)
132 for(
size_t pt=0; pt<quad.size(); pt++) {
134 const Dune::FieldVector<D ,dim>& pos = quad[pt].position();
135 std::vector<Dune::FieldVector<R,1> > q1Values(size);
138 D weight = quad[pt].weight();
140 for (
int k=0; k<size; k++) {
141 integral[k] += q1Values[k]*weight;
143 for (
int l=0; l<=k; l++)
144 massMat[k][l] += weight*(q1Values[k]*q1Values[l]);
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];
155 for (
int i=0; i<size; i++) {
157 Dune::FieldVector<R, size> rhs(0);
158 rhs[i] = integral[i];
161 massMat.solve(coeffs[i] ,rhs);
165 basis.setCoefficients(coeffs);
166 interpolation.setCoefficients(coeffs);
169 template<
class D,
class R,
int dim,
bool faceDual>
170 void DualQ1LocalFiniteElement<D,R,dim,faceDual>::setupFaceDualCoefficients()
173 const int size = 1 <<dim;
174 std::array<Dune::FieldVector<R, size>, size> coeffs;
179 typedef Dune::ReferenceElement<D,dim> RefElement;
180 const RefElement& refElement = Dune::ReferenceElements<D,dim>::general(gt);
183 for (
size_t i=0; i<refElement.size(1);i++) {
185 const Dune::QuadratureRule<D,dim-1>& quad = Dune::QuadratureRules<D,dim-1>::rule(refElement.type(i,1),2*dim);
190 Dune::FieldMatrix<R, size/2, size/2> massMat;
194 const typename RefElement::template Codim<1>::Geometry& geometry = refElement.template geometry<1>(i);
197 std::vector<Dune::FieldVector<R,1> > integral(size/2);
198 for (
int k=0; k<size/2; k++)
201 for(
size_t pt=0; pt<quad.size(); pt++) {
203 const Dune::FieldVector<D,dim-1>& pos = quad[pt].position();
204 const Dune::FieldVector<D,dim> elementPos = geometry.global(pos);
206 std::vector<Dune::FieldVector<R,1> > q1Values(size);
209 D weight = quad[pt].weight();
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;
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]);
225 for (
int l=0; l<refElement.size(i,1,dim); l++) {
227 int row = refElement.subEntity(i,1,l,dim);
228 Dune::FieldVector<R, size/2> rhs(0);
229 rhs[l] = integral[l];
231 Dune::FieldVector<R, size/2> x(0);
232 massMat.solve(x ,rhs);
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];
241 basis.setCoefficients(coeffs);
242 interpolation.setCoefficients(coeffs);
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
Definition: dualq1localinterpolation.hh:17
Definition: brezzidouglasmarini1cube2dlocalbasis.hh:14
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: q1localbasis.hh:37
The local dual Q1 finite element on cubes.
Definition: dualq1.hh:39
Lagrange shape functions of order 1 on the reference cube.
Definition: q1localbasis.hh:24
DualQ1LocalFiniteElement * clone() const
Definition: dualq1.hh:93
traits helper struct
Definition: localfiniteelementtraits.hh:10
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: dualq1.hh:68
LI LocalInterpolationType
Definition: localfiniteelementtraits.hh:22
Layout map for dual Q1 elements.
Definition: dualq1localcoefficients.hh:22
LC LocalCoefficientsType
Definition: localfiniteelementtraits.hh:18
DualQ1LocalFiniteElement()
Definition: dualq1.hh:49
GeometryType type() const
Definition: dualq1.hh:88
LB LocalBasisType
Definition: localfiniteelementtraits.hh:14
const Traits::LocalBasisType & localBasis() const
Definition: dualq1.hh:61
const Traits::LocalInterpolationType & localInterpolation() const
Definition: dualq1.hh:75
Dual Lagrange shape functions of order 1 on the reference cube.
Definition: dualq1localbasis.hh:25