5#ifndef DUNE_GEOMETRY_TEST_LOCALFINITEELEMENT_HH
6#define DUNE_GEOMETRY_TEST_LOCALFINITEELEMENT_HH
17template <
class D,
class R,
unsigned int dim>
18struct ScalarLocalBasisTraits
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>;
33template <
class D,
class R,
unsigned int dim>
37 using Traits = ScalarLocalBasisTraits<D,R,dim>;
40 static constexpr unsigned int size () {
return dim+1; }
43 void evaluateFunction (
const typename Traits::DomainType& x,
44 std::vector<typename Traits::RangeType>& out)
const
48 for (
unsigned int i=0; i<dim; i++)
56 void evaluateJacobian (
const typename Traits::DomainType& x,
57 std::vector<typename Traits::JacobianType>& out)
const
60 std::fill(out[0][0].begin(), out[0][0].end(), -1);
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);
68 static constexpr unsigned int order () {
return 1; }
72template <
class D,
class R,
unsigned int dim>
76 using Traits = ScalarLocalBasisTraits<D,R,dim>;
79 static constexpr unsigned int size () {
return power(2, dim); }
82 void evaluateFunction (
const typename Traits::DomainType& x,
83 std::vector<typename Traits::RangeType>& out)
const
86 for (
size_t i=0; i<
size(); i++)
90 for (
unsigned int j=0; j<dim; j++)
92 out[i] *= (i & (1<<j)) ? x[j] : 1-x[j];
97 void evaluateJacobian (
const typename Traits::DomainType& x,
98 std::vector<typename Traits::JacobianType>& out)
const
103 for (
unsigned int i=0; i<
size(); i++)
106 for (
unsigned int j=0; j<dim; j++)
110 out[i][0][j] = (i & (1<<j)) ? 1 : -1;
112 for (
unsigned int l=0; l<dim; l++)
116 out[i][0][j] *= (i & (1<<l)) ? x[l] : 1-x[l];
123 static constexpr unsigned int order () {
return 1; }
128class P1LocalInterpolation
132 template <
class F,
class C>
135 constexpr auto dim = LB::Traits::dimDomain;
136 out.resize(LB::size());
139 typename LB::Traits::DomainType x;
140 std::fill(x.begin(), x.end(), 0);
144 for (
int i=0; i<dim; i++) {
145 for (
int j=0; j<dim; j++)
154class Q1LocalInterpolation
158 template <
class F,
class C>
161 constexpr auto dim = LB::Traits::dimDomain;
162 out.resize(LB::size());
164 typename LB::Traits::DomainType x;
165 for (
unsigned int i=0; i<LB::size(); i++)
168 for (
int j=0; j<dim; j++)
169 x[j] = (i & (1<<j)) ? 1.0 : 0.0;
178template <
class LB,
template <
class>
class LI>
179class LocalFiniteElement
184 using LocalBasisType = LB;
185 using LocalInterpolationType = LI<LB>;
188 const auto& localBasis ()
const {
return basis_; }
189 const auto& localInterpolation ()
const {
return interpolation_; }
192 static constexpr std::size_t
size () {
return LB::size(); }
196 LI<LB> interpolation_;
199template <
class D,
class R,
int d>
200using P1LocalFiniteElement = LocalFiniteElement<P1LocalBasis<D,R,d>, P1LocalInterpolation>;
202template <
class D,
class R,
int d>
203using Q1LocalFiniteElement = LocalFiniteElement<Q1LocalBasis<D,R,d>, Q1LocalInterpolation>;
207template <
class LFE,
int cdim,
208 class R =
typename LFE::Traits::LocalBasisType::Traits::RangeFieldType>
209class LocalFiniteElementFunction
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>;
219 static_assert(LocalBasis::Traits::dimRange == 1);
222 LocalFiniteElementFunction () =
default;
223 LocalFiniteElementFunction (
const LocalFiniteElement& lfe, std::vector<Range> coefficients)
225 , coefficients_(
std::move(coefficients))
228 Range operator() (
const Domain& local)
const
230 thread_local std::vector<LocalBasisRange> shapeValues;
231 lfe_.localBasis().evaluateFunction(local, shapeValues);
232 assert(shapeValues.size() == coefficients_.size());
234 for (std::size_t i = 0; i < shapeValues.size(); ++i)
235 range.axpy(shapeValues[i], coefficients_[i]);
239 friend auto derivative (
const LocalFiniteElementFunction& f)
241 return [&lfe=f.lfe_,coefficients=f.coefficients_](
const Domain& local) -> Jacobian
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) {
249 shapeJacobians[i].umtv(coefficients[i][j], jacobian[j]);
257 LocalFiniteElement lfe_{};
258 std::vector<Range> coefficients_{};
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