Dune Core Modules (2.8.0)

dualq1.hh
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
10
11#include <dune/geometry/type.hh>
12#include <dune/geometry/referenceelements.hh>
14
15#include <dune/localfunctions/common/localfiniteelementtraits.hh>
16#include <dune/localfunctions/lagrange/lagrangecube.hh>
17#include "dualq1/dualq1localbasis.hh"
18#include "dualq1/dualq1localcoefficients.hh"
19#include "dualq1/dualq1localinterpolation.hh"
20
21namespace Dune
22{
38 template<class D, class R, int dim, bool faceDual=false>
40 {
41 public:
46
50 {
51 if (faceDual)
52 setupFaceDualCoefficients();
53 else
54 setupDualCoefficients();
55 }
56
59 const typename Traits::LocalBasisType& localBasis () const
60 {
61 return basis;
62 }
63
67 {
68 return coefficients;
69 }
70
74 {
75 return interpolation;
76 }
77
79 unsigned int size () const
80 {
81 return basis.size();
82 }
83
86 static constexpr GeometryType type ()
87 {
88 return GeometryTypes::cube(dim);
89 }
90
91 private:
93 void setupFaceDualCoefficients();
94
96 void setupDualCoefficients();
97
101 };
102
103 template<class D, class R, int dim, bool faceDual>
104 void DualQ1LocalFiniteElement<D,R,dim,faceDual>::setupDualCoefficients()
105 {
106
107 const int size = 1 <<dim;
108 std::array<Dune::FieldVector<R, size>, size> coeffs;
109
110 // dual basis functions are linear combinations of Lagrange elements
111 // compute these coefficients here because the basis and the local interpolation needs them
112 const auto& quad = Dune::QuadratureRules<D,dim>::rule(type(), 2*dim);
113
114 // assemble mass matrix on the reference element
116 massMat = 0;
117
118 // and the integrals of the lagrange shape functions
119 std::vector<Dune::FieldVector<R,1> > integral(size);
120 for (int i=0; i<size; i++)
121 integral[i] = 0;
122
123 Dune::Impl::LagrangeCubeLocalBasis<D,R,dim,1> q1Basis;
124 for(size_t pt=0; pt<quad.size(); pt++) {
125
126 const Dune::FieldVector<D ,dim>& pos = quad[pt].position();
127 std::vector<Dune::FieldVector<R,1> > q1Values(size);
128 q1Basis.evaluateFunction(pos,q1Values);
129
130 D weight = quad[pt].weight();
131
132 for (int k=0; k<size; k++) {
133 integral[k] += q1Values[k]*weight;
134
135 for (int l=0; l<=k; l++)
136 massMat[k][l] += weight*(q1Values[k]*q1Values[l]);
137 }
138 }
139
140 // make matrix symmetric
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];
144
145 //solve for the coefficients
146
147 for (int i=0; i<size; i++) {
148
150 rhs[i] = integral[i];
151
152 coeffs[i] = 0;
153 massMat.solve(coeffs[i] ,rhs);
154
155 }
156
157 basis.setCoefficients(coeffs);
158 interpolation.setCoefficients(coeffs);
159 }
160
161 template<class D, class R, int dim, bool faceDual>
162 void DualQ1LocalFiniteElement<D,R,dim,faceDual>::setupFaceDualCoefficients()
163 {
164
165 const int size = 1 <<dim;
166 std::array<Dune::FieldVector<R, size>, size> coeffs;
167
168 // dual basis functions are linear combinations of Lagrange elements
169 Dune::Impl::LagrangeCubeLocalBasis<D,R,dim,1> q1Basis;
170
171 const auto& refElement = Dune::ReferenceElements<D,dim>::general(type());
172
173 // loop over faces
174 for (int i=0; i<refElement.size(1);i++) {
175
176 const auto& quad = Dune::QuadratureRules<D,dim-1>::rule(refElement.type(i,1),2*dim);
177
178 // for each face assemble the mass matrix over the face of all
179 // non-vanishing basis functions,
180 // for cubes refElement.size(i,1,dim)=size/2
181 Dune::FieldMatrix<R, size/2, size/2> massMat;
182 massMat = 0;
183
184 // get geometry
185 const auto& geometry = refElement.template geometry<1>(i);
186
187 // and the integrals of the lagrange shape functions
188 std::vector<Dune::FieldVector<R,1> > integral(size/2);
189 for (int k=0; k<size/2; k++)
190 integral[k] = 0;
191
192 for(size_t pt=0; pt<quad.size(); pt++) {
193
194 const auto& pos = quad[pt].position();
195 const auto& elementPos = geometry.global(pos);
196
197 std::vector<Dune::FieldVector<R,1> > q1Values(size);
198 q1Basis.evaluateFunction(elementPos,q1Values);
199
200 D weight = quad[pt].weight();
201
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;
205
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]);
209 }
210 }
211 }
212
213 // solve for the coefficients
214 // note that we possibly overwrite coefficients for neighbouring faces
215 // this is okay since the coefficients are symmetric
216 for (int l=0; l<refElement.size(i,1,dim); l++) {
217
218 int row = refElement.subEntity(i,1,l,dim);
219 Dune::FieldVector<R, size/2> rhs(0);
220 rhs[l] = integral[l];
221
222 Dune::FieldVector<R, size/2> x(0);
223 massMat.solve(x ,rhs);
224
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];
228 }
229 }
230 }
231
232 basis.setCoefficients(coeffs);
233 interpolation.setCoefficients(coeffs);
234 }
235}
236#endif
void solve(V1 &x, const V2 &b, bool doPivoting=true) const
Solve system A x = b.
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
The local dual Q1 finite element on cubes.
Definition: dualq1.hh:40
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
const Traits::LocalInterpolationType & localInterpolation() const
Definition: dualq1.hh:73
const Traits::LocalBasisType & localBasis() const
Definition: dualq1.hh:59
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: dualq1.hh:66
DualQ1LocalFiniteElement()
Definition: dualq1.hh:49
static constexpr GeometryType type()
Definition: dualq1.hh:86
Definition: dualq1localinterpolation.hh:19
A dense n x m matrix.
Definition: fmatrix.hh:69
vector space out of a tensor product of fields.
Definition: fvector.hh:95
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:123
static const QuadratureRule & rule(const GeometryType &t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
select the appropriate QuadratureRule for GeometryType t and order p
Definition: quadraturerules.hh:280
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
constexpr GeometryType cube(unsigned int dim)
Returns a GeometryType representing a hypercube of dimension dim.
Definition: type.hh:470
Dune namespace.
Definition: alignedallocator.hh:11
static const ReferenceElement & general(const GeometryType &type)
get general reference elements
Definition: referenceelements.hh:196
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
A unique label for each type of element that can occur in a grid.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)