5#ifndef DUNE_DUAL_Q1_LOCALFINITEELEMENT_HH
6#define DUNE_DUAL_Q1_LOCALFINITEELEMENT_HH
10#include <dune/common/fvector.hh>
11#include <dune/common/fmatrix.hh>
13#include <dune/geometry/type.hh>
14#include <dune/geometry/referenceelements.hh>
15#include <dune/geometry/quadraturerules.hh>
40 template<
class D,
class R,
int dim,
bool faceDual=false>
54 setupFaceDualCoefficients();
56 setupDualCoefficients();
88 static constexpr GeometryType
type ()
90 return GeometryTypes::cube(dim);
95 void setupFaceDualCoefficients();
98 void setupDualCoefficients();
105 template<
class D,
class R,
int dim,
bool faceDual>
106 void DualQ1LocalFiniteElement<D,R,dim,faceDual>::setupDualCoefficients()
109 const int size = 1 <<dim;
110 std::array<Dune::FieldVector<R, size>, size> coeffs;
114 const auto& quad = Dune::QuadratureRules<D,dim>::rule(type(), 2*dim);
117 Dune::FieldMatrix<R, size, size> massMat;
121 std::vector<Dune::FieldVector<R,1> > integral(size);
122 for (
int i=0; i<size; i++)
125 Dune::Impl::LagrangeCubeLocalBasis<D,R,dim,1> q1Basis;
126 for(
size_t pt=0; pt<quad.size(); pt++) {
128 const Dune::FieldVector<D ,dim>& pos = quad[pt].position();
129 std::vector<Dune::FieldVector<R,1> > q1Values(size);
130 q1Basis.evaluateFunction(pos,q1Values);
132 D weight = quad[pt].weight();
134 for (
int k=0; k<size; k++) {
135 integral[k] += q1Values[k]*weight;
137 for (
int l=0; l<=k; l++)
138 massMat[k][l] += weight*(q1Values[k]*q1Values[l]);
143 for (
int i=0; i<size-1; i++)
144 for (
int j=i+1; j<size; j++)
145 massMat[i][j] = massMat[j][i];
149 for (
int i=0; i<size; i++) {
151 Dune::FieldVector<R, size> rhs(0);
152 rhs[i] = integral[i];
155 massMat.solve(coeffs[i] ,rhs);
159 basis.setCoefficients(coeffs);
160 interpolation.setCoefficients(coeffs);
163 template<
class D,
class R,
int dim,
bool faceDual>
164 void DualQ1LocalFiniteElement<D,R,dim,faceDual>::setupFaceDualCoefficients()
167 const int size = 1 <<dim;
168 std::array<Dune::FieldVector<R, size>, size> coeffs;
171 Dune::Impl::LagrangeCubeLocalBasis<D,R,dim,1> q1Basis;
173 const auto& refElement = Dune::ReferenceElements<D,dim>::general(type());
176 for (
int i=0; i<refElement.size(1);i++) {
178 const auto& quad = Dune::QuadratureRules<D,dim-1>::rule(refElement.type(i,1),2*dim);
183 Dune::FieldMatrix<R, size/2, size/2> massMat;
187 const auto& geometry = refElement.template geometry<1>(i);
190 std::vector<Dune::FieldVector<R,1> > integral(size/2);
191 for (
int k=0; k<size/2; k++)
194 for(
size_t pt=0; pt<quad.size(); pt++) {
196 const auto& pos = quad[pt].position();
197 const auto& elementPos = geometry.global(pos);
199 std::vector<Dune::FieldVector<R,1> > q1Values(size);
200 q1Basis.evaluateFunction(elementPos,q1Values);
202 D weight = quad[pt].weight();
204 for (
int k=0; k<refElement.size(i,1,dim); k++) {
205 int row = refElement.subEntity(i,1,k,dim);
206 integral[k] += q1Values[row]*weight;
208 for (
int l=0; l<refElement.size(i,1,dim); l++) {
209 int col = refElement.subEntity(i,1,l,dim);
210 massMat[k][l] += weight*(q1Values[row]*q1Values[col]);
218 for (
int l=0; l<refElement.size(i,1,dim); l++) {
220 int row = refElement.subEntity(i,1,l,dim);
221 Dune::FieldVector<R, size/2> rhs(0);
222 rhs[l] = integral[l];
224 Dune::FieldVector<R, size/2> x(0);
225 massMat.solve(x ,rhs);
227 for (
int k=0; k<refElement.size(i,1,dim); k++) {
228 int col = refElement.subEntity(i,1,k,dim);
229 coeffs[row][col]=x[k];
234 basis.setCoefficients(coeffs);
235 interpolation.setCoefficients(coeffs);
Definition bdfmcube.hh:18
traits helper struct
Definition localfiniteelementtraits.hh:13
LB LocalBasisType
Definition localfiniteelementtraits.hh:16
LC LocalCoefficientsType
Definition localfiniteelementtraits.hh:20
LI LocalInterpolationType
Definition localfiniteelementtraits.hh:24
The local dual Q1 finite element on cubes.
Definition dualq1.hh:42
LocalFiniteElementTraits< DualQ1LocalBasis< D, R, dim >, DualQ1LocalCoefficients< dim >, DualQ1LocalInterpolation< dim, DualQ1LocalBasis< D, R, dim > > > Traits
Definition dualq1.hh:47
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:51
static constexpr GeometryType type()
Definition dualq1.hh:88
Dual Lagrange shape functions of order 1 on the reference cube.
Definition dualq1localbasis.hh:30
Layout map for dual Q1 elements.
Definition dualq1localcoefficients.hh:26
Definition dualq1localinterpolation.hh:22