DUNE-FEM (unstable)

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// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
4// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5#ifndef DUNE_DUAL_Q1_LOCALFINITEELEMENT_HH
6#define DUNE_DUAL_Q1_LOCALFINITEELEMENT_HH
7
8#include <array>
9
12
13#include <dune/geometry/type.hh>
14#include <dune/geometry/referenceelements.hh>
16
17#include <dune/localfunctions/common/localfiniteelementtraits.hh>
18#include <dune/localfunctions/lagrange/lagrangecube.hh>
19#include "dualq1/dualq1localbasis.hh"
20#include "dualq1/dualq1localcoefficients.hh"
21#include "dualq1/dualq1localinterpolation.hh"
22
23namespace Dune
24{
40 template<class D, class R, int dim, bool faceDual=false>
42 {
43 public:
48
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 static constexpr GeometryType type ()
89 {
90 return GeometryTypes::cube(dim);
91 }
92
93 private:
95 void setupFaceDualCoefficients();
96
98 void setupDualCoefficients();
99
101 DualQ1LocalCoefficients<dim> coefficients;
103 };
104
105 template<class D, class R, int dim, bool faceDual>
106 void DualQ1LocalFiniteElement<D,R,dim,faceDual>::setupDualCoefficients()
107 {
108
109 const int size = 1 <<dim;
110 std::array<Dune::FieldVector<R, size>, size> coeffs;
111
112 // dual basis functions are linear combinations of Lagrange elements
113 // compute these coefficients here because the basis and the local interpolation needs them
114 const auto& quad = Dune::QuadratureRules<D,dim>::rule(type(), 2*dim);
115
116 // assemble mass matrix on the reference element
118 massMat = 0;
119
120 // and the integrals of the lagrange shape functions
121 std::vector<Dune::FieldVector<R,1> > integral(size);
122 for (int i=0; i<size; i++)
123 integral[i] = 0;
124
125 Dune::Impl::LagrangeCubeLocalBasis<D,R,dim,1> q1Basis;
126 for(size_t pt=0; pt<quad.size(); pt++) {
127
128 const Dune::FieldVector<D ,dim>& pos = quad[pt].position();
129 std::vector<Dune::FieldVector<R,1> > q1Values(size);
130 q1Basis.evaluateFunction(pos,q1Values);
131
132 D weight = quad[pt].weight();
133
134 for (int k=0; k<size; k++) {
135 integral[k] += q1Values[k]*weight;
136
137 for (int l=0; l<=k; l++)
138 massMat[k][l] += weight*(q1Values[k]*q1Values[l]);
139 }
140 }
141
142 // make matrix symmetric
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];
146
147 //solve for the coefficients
148
149 for (int i=0; i<size; i++) {
150
152 rhs[i] = integral[i];
153
154 coeffs[i] = 0;
155 massMat.solve(coeffs[i] ,rhs);
156
157 }
158
159 basis.setCoefficients(coeffs);
160 interpolation.setCoefficients(coeffs);
161 }
162
163 template<class D, class R, int dim, bool faceDual>
164 void DualQ1LocalFiniteElement<D,R,dim,faceDual>::setupFaceDualCoefficients()
165 {
166
167 const int size = 1 <<dim;
168 std::array<Dune::FieldVector<R, size>, size> coeffs;
169
170 // dual basis functions are linear combinations of Lagrange elements
171 Dune::Impl::LagrangeCubeLocalBasis<D,R,dim,1> q1Basis;
172
173 const auto& refElement = Dune::ReferenceElements<D,dim>::general(type());
174
175 // loop over faces
176 for (int i=0; i<refElement.size(1);i++) {
177
178 const auto& quad = Dune::QuadratureRules<D,dim-1>::rule(refElement.type(i,1),2*dim);
179
180 // for each face assemble the mass matrix over the face of all
181 // non-vanishing basis functions,
182 // for cubes refElement.size(i,1,dim)=size/2
183 Dune::FieldMatrix<R, size/2, size/2> massMat;
184 massMat = 0;
185
186 // get geometry
187 const auto& geometry = refElement.template geometry<1>(i);
188
189 // and the integrals of the lagrange shape functions
190 std::vector<Dune::FieldVector<R,1> > integral(size/2);
191 for (int k=0; k<size/2; k++)
192 integral[k] = 0;
193
194 for(size_t pt=0; pt<quad.size(); pt++) {
195
196 const auto& pos = quad[pt].position();
197 const auto& elementPos = geometry.global(pos);
198
199 std::vector<Dune::FieldVector<R,1> > q1Values(size);
200 q1Basis.evaluateFunction(elementPos,q1Values);
201
202 D weight = quad[pt].weight();
203
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;
207
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]);
211 }
212 }
213 }
214
215 // solve for the coefficients
216 // note that we possibly overwrite coefficients for neighbouring faces
217 // this is okay since the coefficients are symmetric
218 for (int l=0; l<refElement.size(i,1,dim); l++) {
219
220 int row = refElement.subEntity(i,1,l,dim);
221 Dune::FieldVector<R, size/2> rhs(0);
222 rhs[l] = integral[l];
223
224 Dune::FieldVector<R, size/2> x(0);
225 massMat.solve(x ,rhs);
226
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];
230 }
231 }
232 }
233
234 basis.setCoefficients(coeffs);
235 interpolation.setCoefficients(coeffs);
236 }
237}
238#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:30
Layout map for dual Q1 elements.
Definition: dualq1localcoefficients.hh:26
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
Definition: dualq1localinterpolation.hh:22
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:114
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:326
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:462
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
static const ReferenceElement & general(const GeometryType &type)
get general reference elements
Definition: referenceelements.hh:156
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
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 12, 23:30, 2024)