DUNE-FEM (unstable)

localfiniteelements.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_GEOMETRY_TEST_LOCALFINITEELEMENT_HH
6#define DUNE_GEOMETRY_TEST_LOCALFINITEELEMENT_HH
7
8#include <algorithm>
9#include <vector>
10
13#include <dune/common/math.hh>
14
15namespace Dune::Impl {
16
17template <class D, class R, unsigned int dim>
18struct ScalarLocalBasisTraits
19{
20 using DomainFieldType = D;
21 using RangeFieldType = R;
22 using DomainType = FieldVector<D,dim>;
23 using RangeType = FieldVector<R,1>;
24 using JacobianType = FieldMatrix<R,1,dim>;
25
26 enum {
27 dimDomain = dim,
28 dimRange = 1
29 };
30};
31
33template <class D, class R, unsigned int dim>
34class P1LocalBasis
35{
36public:
37 using Traits = ScalarLocalBasisTraits<D,R,dim>;
38
40 static constexpr unsigned int size () { return dim+1; }
41
43 void evaluateFunction (const typename Traits::DomainType& x,
44 std::vector<typename Traits::RangeType>& out) const
45 {
46 out.resize(size());
47 out[0] = 1.0;
48 for (unsigned int i=0; i<dim; i++)
49 {
50 out[0] -= x[i];
51 out[i+1] = x[i];
52 }
53 }
54
56 void evaluateJacobian (const typename Traits::DomainType& x,
57 std::vector<typename Traits::JacobianType>& out) const
58 {
59 out.resize(size());
60 std::fill(out[0][0].begin(), out[0][0].end(), -1);
61
62 for (unsigned int i=0; i<dim; i++)
63 for (unsigned int j=0; j<dim; j++)
64 out[i+1][0][j] = (i==j);
65 }
66
68 static constexpr unsigned int order () { return 1; }
69};
70
72template <class D, class R, unsigned int dim>
73class Q1LocalBasis
74{
75public:
76 using Traits = ScalarLocalBasisTraits<D,R,dim>;
77
79 static constexpr unsigned int size () { return power(2, dim); }
80
82 void evaluateFunction (const typename Traits::DomainType& x,
83 std::vector<typename Traits::RangeType>& out) const
84 {
85 out.resize(size());
86 for (size_t i=0; i<size(); i++)
87 {
88 out[i] = 1;
89
90 for (unsigned int j=0; j<dim; j++)
91 // if j-th bit of i is set multiply with x[j], else with 1-x[j]
92 out[i] *= (i & (1<<j)) ? x[j] : 1-x[j];
93 }
94 }
95
97 void evaluateJacobian (const typename Traits::DomainType& x,
98 std::vector<typename Traits::JacobianType>& out) const
99 {
100 out.resize(size());
101
102 // Loop over all shape functions
103 for (unsigned int i=0; i<size(); i++)
104 {
105 // Loop over all coordinate directions
106 for (unsigned int j=0; j<dim; j++)
107 {
108 // Initialize: the overall expression is a product
109 // if j-th bit of i is set to 1, else -11
110 out[i][0][j] = (i & (1<<j)) ? 1 : -1;
111
112 for (unsigned int l=0; l<dim; l++)
113 {
114 if (j!=l)
115 // if l-th bit of i is set multiply with x[l], else with 1-x[l]
116 out[i][0][j] *= (i & (1<<l)) ? x[l] : 1-x[l];
117 }
118 }
119 }
120 }
121
123 static constexpr unsigned int order () { return 1; }
124};
125
126
127template <class LB>
128class P1LocalInterpolation
129{
130public:
132 template <class F, class C>
133 void interpolate (F f, std::vector<C>& out) const
134 {
135 constexpr auto dim = LB::Traits::dimDomain;
136 out.resize(LB::size());
137
138 // vertex 0
139 typename LB::Traits::DomainType x;
140 std::fill(x.begin(), x.end(), 0);
141 out[0] = f(x);
142
143 // remaining vertices
144 for (int i=0; i<dim; i++) {
145 for (int j=0; j<dim; j++)
146 x[j] = (i==j);
147
148 out[i+1] = f(x);
149 }
150 }
151};
152
153template <class LB>
154class Q1LocalInterpolation
155{
156public:
158 template <class F, class C>
159 void interpolate (F f, std::vector<C>& out) const
160 {
161 constexpr auto dim = LB::Traits::dimDomain;
162 out.resize(LB::size());
163
164 typename LB::Traits::DomainType x;
165 for (unsigned int i=0; i<LB::size(); i++)
166 {
167 // Generate coordinate of the i-th corner of the reference cube
168 for (int j=0; j<dim; j++)
169 x[j] = (i & (1<<j)) ? 1.0 : 0.0;
170
171 out[i] = f(x);
172 }
173 }
174};
175
176
178template <class LB, template <class> class LI>
179class LocalFiniteElement
180{
181public:
182 struct Traits
183 {
184 using LocalBasisType = LB;
185 using LocalInterpolationType = LI<LB>;
186 };
187
188 const auto& localBasis () const { return basis_; }
189 const auto& localInterpolation () const { return interpolation_; }
190
192 static constexpr std::size_t size () { return LB::size(); }
193
194private:
195 LB basis_;
196 LI<LB> interpolation_;
197};
198
199template <class D, class R, int d>
200using P1LocalFiniteElement = LocalFiniteElement<P1LocalBasis<D,R,d>, P1LocalInterpolation>;
201
202template <class D, class R, int d>
203using Q1LocalFiniteElement = LocalFiniteElement<Q1LocalBasis<D,R,d>, Q1LocalInterpolation>;
204
205
206
207template <class LFE, int cdim,
208 class R = typename LFE::Traits::LocalBasisType::Traits::RangeFieldType>
209class LocalFiniteElementFunction
210{
211 using LocalFiniteElement = LFE;
212 using LocalBasis = typename LocalFiniteElement::Traits::LocalBasisType;
213 using LocalBasisRange = typename LocalBasis::Traits::RangeType;
214 using LocalBasisJacobian = typename LocalBasis::Traits::JacobianType;
215 using Domain = typename LocalBasis::Traits::DomainType;
216 using Range = FieldVector<R,cdim>;
217 using Jacobian = FieldMatrix<R,cdim,LocalBasis::Traits::dimDomain>;
218
219 static_assert(LocalBasis::Traits::dimRange == 1);
220
221public:
222 LocalFiniteElementFunction () = default;
223 LocalFiniteElementFunction (const LocalFiniteElement& lfe, std::vector<Range> coefficients)
224 : lfe_(lfe)
225 , coefficients_(std::move(coefficients))
226 {}
227
228 Range operator() (const Domain& local) const
229 {
230 thread_local std::vector<LocalBasisRange> shapeValues;
231 lfe_.localBasis().evaluateFunction(local, shapeValues);
232 assert(shapeValues.size() == coefficients_.size());
233 Range range(0);
234 for (std::size_t i = 0; i < shapeValues.size(); ++i)
235 range.axpy(shapeValues[i], coefficients_[i]);
236 return range;
237 }
238
239 friend auto derivative (const LocalFiniteElementFunction& f)
240 {
241 return [&lfe=f.lfe_,coefficients=f.coefficients_](const Domain& local) -> Jacobian
242 {
243 thread_local std::vector<LocalBasisJacobian> shapeJacobians;
244 lfe.localBasis().evaluateJacobian(local, shapeJacobians);
245 assert(shapeJacobians.size() == coefficients.size());
246 Jacobian jacobian(0);
247 for (std::size_t i = 0; i < shapeJacobians.size(); ++i) {
248 for (int j = 0; j < Jacobian::rows; ++j) {
249 shapeJacobians[i].umtv(coefficients[i][j], jacobian[j]);
250 }
251 }
252 return jacobian;
253 };
254 }
255
256private:
257 LocalFiniteElement lfe_{};
258 std::vector<Range> coefficients_{};
259};
260
261} // end namespace Dune::Impl
262
263#endif // DUNE_GEOMETRY_TEST_LOCALFINITEELEMENT_HH
static constexpr int rows
The number of rows.
Definition: fmatrix.hh:123
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.
static void interpolate(const GridFunction &u, DiscreteFunction &v)
perform native interpolation of a discrete function space
Definition: interpolate.hh:54
Some useful basic math stuff.
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
constexpr Base power(Base m, Exponent p)
Power method for integer exponents.
Definition: math.hh:75
STL namespace.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)