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>
52 setupFaceDualCoefficients();
54 setupDualCoefficients();
86 static constexpr GeometryType
type ()
88 return GeometryTypes::cube(dim);
93 void setupFaceDualCoefficients();
96 void setupDualCoefficients();
99 DualQ1LocalCoefficients<dim> coefficients;
103 template<
class D,
class R,
int dim,
bool faceDual>
107 const int size = 1 <<dim;
108 std::array<Dune::FieldVector<R, size>, size> coeffs;
112 const auto& quad = Dune::QuadratureRules<D,dim>::rule(
type(), 2*dim);
115 Dune::FieldMatrix<R, size, size> massMat;
119 std::vector<Dune::FieldVector<R,1> > integral(size);
120 for (
int i=0; i<
size; i++)
124 for(
size_t pt=0; pt<quad.size(); pt++) {
126 const Dune::FieldVector<D ,dim>& pos = quad[pt].position();
127 std::vector<Dune::FieldVector<R,1> > q1Values(size);
130 D weight = quad[pt].weight();
132 for (
int k=0; k<
size; k++) {
133 integral[k] += q1Values[k]*weight;
135 for (
int l=0; l<=k; l++)
136 massMat[k][l] += weight*(q1Values[k]*q1Values[l]);
141 for (
int i=0; i<size-1; i++)
142 for (
int j=i+1; j<
size; j++)
143 massMat[i][j] = massMat[j][i];
147 for (
int i=0; i<
size; i++) {
149 Dune::FieldVector<R, size> rhs(0);
150 rhs[i] = integral[i];
153 massMat.solve(coeffs[i] ,rhs);
157 basis.setCoefficients(coeffs);
158 interpolation.setCoefficients(coeffs);
161 template<
class D,
class R,
int dim,
bool faceDual>
165 const int size = 1 <<dim;
166 std::array<Dune::FieldVector<R, size>, size> coeffs;
171 const auto& refElement = Dune::ReferenceElements<D,dim>::general(
type());
174 for (
size_t i=0; i<refElement.size(1);i++) {
176 const auto& quad = Dune::QuadratureRules<D,dim-1>::rule(refElement.type(i,1),2*dim);
181 Dune::FieldMatrix<R, size/2, size/2> massMat;
185 const auto& geometry = refElement.template geometry<1>(i);
188 std::vector<Dune::FieldVector<R,1> > integral(size/2);
189 for (
int k=0; k<size/2; k++)
192 for(
size_t pt=0; pt<quad.size(); pt++) {
194 const auto& pos = quad[pt].position();
195 const auto& elementPos = geometry.global(pos);
197 std::vector<Dune::FieldVector<R,1> > q1Values(size);
200 D weight = quad[pt].weight();
202 for (
int k=0; k<refElement.size(i,1,dim); k++) {
203 int row = refElement.subEntity(i,1,k,dim);
204 integral[k] += q1Values[row]*weight;
206 for (
int l=0; l<refElement.size(i,1,dim); l++) {
207 int col = refElement.subEntity(i,1,l,dim);
208 massMat[k][l] += weight*(q1Values[row]*q1Values[col]);
216 for (
int l=0; l<refElement.size(i,1,dim); l++) {
218 int row = refElement.subEntity(i,1,l,dim);
219 Dune::FieldVector<R, size/2> rhs(0);
220 rhs[l] = integral[l];
222 Dune::FieldVector<R, size/2> x(0);
223 massMat.solve(x ,rhs);
225 for (
int k=0; k<refElement.size(i,1,dim); k++) {
226 int col = refElement.subEntity(i,1,k,dim);
227 coeffs[row][col]=x[k];
232 basis.setCoefficients(coeffs);
233 interpolation.setCoefficients(coeffs);
traits helper struct
Definition: localfiniteelementtraits.hh:10
static constexpr GeometryType type()
Definition: dualq1.hh:86
The local dual Q1 finite element on cubes.
Definition: dualq1.hh:39
const Traits::LocalInterpolationType & localInterpolation() const
Definition: dualq1.hh:73
LI LocalInterpolationType
Definition: localfiniteelementtraits.hh:22
Dual Lagrange shape functions of order 1 on the reference cube.
Definition: dualq1localbasis.hh:26
LC LocalCoefficientsType
Definition: localfiniteelementtraits.hh:18
Definition: dualq1localinterpolation.hh:17
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: q1localbasis.hh:39
const Traits::LocalBasisType & localBasis() const
Definition: dualq1.hh:59
Definition: brezzidouglasmarini1cube2dlocalbasis.hh:15
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: dualq1.hh:66
DualQ1LocalFiniteElement()
Definition: dualq1.hh:49
LB LocalBasisType
Definition: localfiniteelementtraits.hh:14
Lagrange shape functions of order 1 on the reference cube.
Definition: q1localbasis.hh:26
Layout map for dual Q1 elements.
Definition: dualq1localcoefficients.hh:22
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:79