dune-localfunctions  2.9.0
dualq1.hh
Go to the documentation of this file.
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 (C) 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 
10 #include <dune/common/fvector.hh>
11 #include <dune/common/fmatrix.hh>
12 
13 #include <dune/geometry/type.hh>
14 #include <dune/geometry/referenceelements.hh>
15 #include <dune/geometry/quadraturerules.hh>
16 
22 
23 namespace 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
117  Dune::FieldMatrix<R, size, size> massMat;
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 
151  Dune::FieldVector<R, size> rhs(0);
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
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
Dual Lagrange shape functions of order 1 on the reference cube.
Definition: dualq1localbasis.hh:29
Layout map for dual Q1 elements.
Definition: dualq1localcoefficients.hh:25
Definition: dualq1localinterpolation.hh:21
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
DualQ1LocalFiniteElement()
Definition: dualq1.hh:51
const Traits::LocalBasisType & localBasis() const
Definition: dualq1.hh:61
const Traits::LocalInterpolationType & localInterpolation() const
Definition: dualq1.hh:75
static constexpr GeometryType type()
Definition: dualq1.hh:88
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: dualq1.hh:68