│ │ │ │
│ │ │ │
│ │ │ │
│ │ │ │
│ │ │ │
│ │ │ │ -
7#ifndef DUNE_FUNCTIONS_GRIDFUNCTIONS_GRID_FUNCTION_IMP_HH
│ │ │ │ -
8#define DUNE_FUNCTIONS_GRIDFUNCTIONS_GRID_FUNCTION_IMP_HH
│ │ │ │ +
7#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH
│ │ │ │ +
8#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
18#include <dune/common/dynmatrix.hh>
│ │ │ │
│ │ │ │ -
23struct HasFreeLocalFunction
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
26 auto require(F&& f) ->
decltype(
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
20#include <dune/localfunctions/common/localbasis.hh>
│ │ │ │ +
21#include <dune/common/diagonalmatrix.hh>
│ │ │ │ +
22#include <dune/localfunctions/common/localkey.hh>
│ │ │ │ +
23#include <dune/localfunctions/common/localfiniteelementtraits.hh>
│ │ │ │ +
24#include <dune/geometry/type.hh>
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
37template<
class Signature,
class DerivativeInterface,
class LocalFunctionInterface,
class EntitySet>
│ │ │ │ -
38class GridFunctionWrapperInterface :
│ │ │ │ -
39 public DifferentiableFunctionWrapperInterface<Signature, DerivativeInterface>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
42 virtual LocalFunctionInterface wrappedLocalFunction()
const = 0;
│ │ │ │ -
│ │ │ │ -
44 virtual const EntitySet& wrappedEntitySet()
const = 0;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
49template<
class Signature,
class DerivativeInterface,
class LocalFunctionInterface,
class EntitySet,
class B>
│ │ │ │ -
50class GridFunctionWrapperImplementation :
│ │ │ │ -
51 public DifferentiableFunctionWrapperImplementation<Signature, DerivativeInterface, B>
│ │ │ │ -
│ │ │ │ -
53 using Base = DifferentiableFunctionWrapperImplementation<Signature, DerivativeInterface, B>;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
57 virtual LocalFunctionInterface wrappedLocalFunction()
const
│ │ │ │ -
│ │ │ │ -
59 return localFunction(this->get());
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
34template<
typename GV,
typename R>
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
49template<
class GV,
class R>
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
54 typedef typename GV::ctype D;
│ │ │ │ +
55 enum {dim = GV::dimension};
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
59 typedef LocalBasisTraits<D,dim,FieldVector<D,dim>,R,1,FieldVector<R,1>,
│ │ │ │ +
│ │ │ │
│ │ │ │ -
62 virtual const EntitySet& wrappedEntitySet()
const
│ │ │ │ -
│ │ │ │ -
64 return this->get().entitySet();
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
68 : preBasis_(preBasis),
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
76 std::vector<FieldVector<R,1> >& out)
const
│ │ │ │ +
│ │ │ │ +
78 FieldVector<D,dim> globalIn = offset_;
│ │ │ │ +
79 scaling_.umv(in,globalIn);
│ │ │ │ +
│ │ │ │ +
81 preBasis_.evaluateFunction(globalIn, out, lFE_.currentKnotSpan_);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
88 std::vector<FieldMatrix<D,1,dim> >& out)
const
│ │ │ │ +
│ │ │ │ +
90 FieldVector<D,dim> globalIn = offset_;
│ │ │ │ +
91 scaling_.umv(in,globalIn);
│ │ │ │ +
│ │ │ │ +
93 preBasis_.evaluateJacobian(globalIn, out, lFE_.currentKnotSpan_);
│ │ │ │ +
│ │ │ │ +
95 for (
size_t i=0; i<out.size(); i++)
│ │ │ │ +
96 for (
int j=0; j<dim; j++)
│ │ │ │ +
97 out[i][0][j] *= scaling_[j][j];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
102 inline void evaluate (
const typename std::array<int,k>& directions,
│ │ │ │ +
103 const typename Traits::DomainType& in,
│ │ │ │ +
104 std::vector<typename Traits::RangeType>& out)
const
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
113 FieldVector<D,dim> globalIn = offset_;
│ │ │ │ +
114 scaling_.umv(in,globalIn);
│ │ │ │ +
│ │ │ │ +
116 preBasis_.evaluate(directions, globalIn, out, lFE_.currentKnotSpan_);
│ │ │ │ +
│ │ │ │ +
118 for (
size_t i=0; i<out.size(); i++)
│ │ │ │ +
119 out[i][0] *= scaling_[directions[0]][directions[0]];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
124 FieldVector<D,dim> globalIn = offset_;
│ │ │ │ +
125 scaling_.umv(in,globalIn);
│ │ │ │ +
│ │ │ │ +
127 preBasis_.evaluate(directions, globalIn, out, lFE_.currentKnotSpan_);
│ │ │ │ +
│ │ │ │ +
129 for (
size_t i=0; i<out.size(); i++)
│ │ │ │ +
130 out[i][0] *= scaling_[directions[0]][directions[0]]*scaling_[directions[1]][directions[1]];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
134 DUNE_THROW(NotImplemented,
"B-Spline derivatives of order " << k <<
" not implemented yet!");
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
147 return *std::max_element(preBasis_.order_.begin(), preBasis_.order_.end());
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
164 FieldVector<D,dim> offset_;
│ │ │ │ +
165 DiagonalMatrix<D,dim> scaling_;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
185 std::array<unsigned int,dim> multiindex (
unsigned int i)
const
│ │ │ │ +
│ │ │ │ +
187 std::array<unsigned int,dim> alpha;
│ │ │ │ +
188 for (
int j=0; j<dim; j++)
│ │ │ │ +
│ │ │ │ +
190 alpha[j] = i % sizes_[j];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
197 void setup1d(std::vector<unsigned int>& subEntity)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
208 unsigned lastIndex=0;
│ │ │ │ +
209 subEntity[lastIndex++] = 0;
│ │ │ │ +
210 for (
unsigned i = 0; i < sizes_[0] - 2; ++i)
│ │ │ │ +
211 subEntity[lastIndex++] = 0;
│ │ │ │ +
│ │ │ │ +
213 subEntity[lastIndex++] = 1;
│ │ │ │ +
│ │ │ │ +
215 assert(
size()==lastIndex);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
218 void setup2d(std::vector<unsigned int>& subEntity)
│ │ │ │ +
│ │ │ │ +
220 unsigned lastIndex=0;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
234 subEntity[lastIndex++] = 0;
│ │ │ │ +
235 for (
unsigned i = 0; i < sizes_[0]-2; ++i)
│ │ │ │ +
236 subEntity[lastIndex++] = 2;
│ │ │ │ +
│ │ │ │ +
238 subEntity[lastIndex++] = 1;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
241 for (
unsigned e = 0; e < sizes_[1]-2; ++e)
│ │ │ │ +
│ │ │ │ +
243 subEntity[lastIndex++] = 0;
│ │ │ │ +
244 for (
unsigned i = 0; i < sizes_[0]-2; ++i)
│ │ │ │ +
245 subEntity[lastIndex++] = 0;
│ │ │ │ +
246 subEntity[lastIndex++] = 1;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
250 subEntity[lastIndex++] = 2;
│ │ │ │ +
251 for (
unsigned i = 0; i < sizes_[0]-2; ++i)
│ │ │ │ +
252 subEntity[lastIndex++] = 3;
│ │ │ │ +
│ │ │ │ +
254 subEntity[lastIndex++] = 3;
│ │ │ │ +
│ │ │ │ +
256 assert(
size()==lastIndex);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
261 void init(
const std::array<unsigned,dim>& sizes)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
268 std::vector<unsigned int> codim(li_.size());
│ │ │ │ +
│ │ │ │ +
270 for (std::size_t i=0; i<codim.size(); i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
275 std::array<unsigned int,dim> mIdx = multiindex(i);
│ │ │ │ +
276 for (
int j=0; j<dim; j++)
│ │ │ │ +
277 if (mIdx[j]==0 or mIdx[j]==sizes[j]-1)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
286 std::vector<unsigned int> index(
size());
│ │ │ │ +
│ │ │ │ +
288 for (std::size_t i=0; i<index.size(); i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
292 std::array<unsigned int,dim> mIdx = multiindex(i);
│ │ │ │ +
│ │ │ │ +
294 for (
int j=dim-1; j>=0; j--)
│ │ │ │ +
295 if (mIdx[j]>0 and mIdx[j]<sizes[j]-1)
│ │ │ │ +
296 index[i] = (sizes[j]-1)*index[i] + (mIdx[j]-1);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
300 std::vector<unsigned int> subEntity(li_.size());
│ │ │ │ +
│ │ │ │ +
302 if (subEntity.size() > 0)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
308 }
else if (dim==2 and sizes_[0]>1 and sizes_[1]>1) {
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
315 for (
size_t i=0; i<li_.size(); i++)
│ │ │ │ +
316 li_[i] = LocalKey(subEntity[i], codim[i], index[i]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
322 return std::accumulate(sizes_.begin(), sizes_.end(), 1, std::multiplies<unsigned int>());
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
334 std::array<unsigned, dim> sizes_;
│ │ │ │ +
│ │ │ │ +
336 std::vector<LocalKey> li_;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
343template<
int dim,
class LB>
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
348 template<
typename F,
typename C>
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
351 DUNE_THROW(NotImplemented,
"BSplineLocalInterpolation::interpolate");
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
365template<
class GV,
class R>
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
368 typedef typename GV::ctype D;
│ │ │ │ +
369 enum {dim = GV::dimension};
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
375 typedef LocalFiniteElementTraits<BSplineLocalBasis<GV,R>,
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
399 void bind(
const std::array<unsigned,dim>& elementIdx)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
402 for (
size_t i=0; i<elementIdx.size(); i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
410 for (
size_t j=0; j<elementIdx[i]; j++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
425 std::array<unsigned int, dim> sizes;
│ │ │ │ +
426 for (
size_t i=0; i<dim; i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
453 for (
int i=0; i<dim; i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
462 return GeometryTypes::cube(dim);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
471 unsigned int r = order[i]+1;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
508 static const int dim = GV::dimension;
│ │ │ │ +
│ │ │ │ +
511 class MultiDigitCounter
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
518 MultiDigitCounter(
const std::array<unsigned int,dim>& limits)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
521 std::fill(counter_.begin(), counter_.end(), 0);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
525 MultiDigitCounter& operator++()
│ │ │ │ +
│ │ │ │ +
527 for (
int i=0; i<dim; i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
532 if (counter_[i] < limits_[i])
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
541 const unsigned int& operator[](
int i)
const
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
547 unsigned int cycle()
const
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
550 for (
int i=0; i<dim; i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
558 const std::array<unsigned int,dim> limits_;
│ │ │ │ +
│ │ │ │ +
561 std::array<unsigned int,dim> counter_;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
595 const std::vector<double>& knotVector,
│ │ │ │ +
│ │ │ │ +
597 bool makeOpen =
true)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
607 for (
int i=0; i<dim; i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
612 for (
unsigned int j=0; j<order; j++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
618 for (
unsigned int j=0; j<order; j++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
647 const FieldVector<double,dim>& lowerLeft,
│ │ │ │ +
648 const FieldVector<double,dim>& upperRight,
│ │ │ │ +
649 const std::array<unsigned int,dim>& elements,
│ │ │ │ +
│ │ │ │ +
651 bool makeOpen =
true)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
659 for (
int i=0; i<dim; i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
664 for (
unsigned int j=0; j<order; j++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
668 for (
size_t j=0; j<elements[i]+1; j++)
│ │ │ │ +
669 knotVectors_[i].push_back(lowerLeft[i] + j*(upperRight[i]-lowerLeft[i]) / elements[i]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
672 for (
unsigned int j=0; j<order; j++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
707 for (
int i=0; i<dim; i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
713 template<
typename It>
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
718 std::array<unsigned int, dim> localSizes;
│ │ │ │ +
719 for (
int i=0; i<dim; i++)
│ │ │ │ +
│ │ │ │ +
721 for (
size_type i = 0, end = node.
size() ; i < end ; ++i, ++it)
│ │ │ │ +
│ │ │ │ +
723 std::array<unsigned int,dim> localIJK =
getIJK(i, localSizes);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
726 const auto order =
order_;
│ │ │ │ +
│ │ │ │ +
728 std::array<unsigned int,dim> globalIJK;
│ │ │ │ +
729 for (
int i=0; i<dim; i++)
│ │ │ │ +
730 globalIJK[i] = std::max((
int)currentKnotSpan[i] - (
int)order[i], 0) + localIJK[i];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
735 for (
int i=dim-2; i>=0; i--)
│ │ │ │ +
736 globalIdx = globalIdx *
size(i) + globalIJK[i];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
746 unsigned int result = 1;
│ │ │ │ +
747 for (
size_t i=0; i<dim; i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
753 unsigned int size (
size_t d)
const
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
763 std::vector<FieldVector<R,1> >& out,
│ │ │ │ +
764 const std::array<unsigned,dim>& currentKnotSpan)
const
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
767 std::array<std::vector<R>, dim> oneDValues;
│ │ │ │ +
│ │ │ │ +
769 for (
size_t i=0; i<dim; i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
772 std::array<unsigned int, dim> limits;
│ │ │ │ +
773 for (
int i=0; i<dim; i++)
│ │ │ │ +
774 limits[i] = oneDValues[i].
size();
│ │ │ │ +
│ │ │ │ +
776 MultiDigitCounter ijkCounter(limits);
│ │ │ │ +
│ │ │ │ +
778 out.resize(ijkCounter.cycle());
│ │ │ │ +
│ │ │ │ +
780 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
783 for (
size_t j=0; j<dim; j++)
│ │ │ │ +
784 out[i] *= oneDValues[j][ijkCounter[j]];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
794 std::vector<FieldMatrix<R,1,dim> >& out,
│ │ │ │ +
795 const std::array<unsigned,dim>& currentKnotSpan)
const
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
798 std::array<unsigned int, dim> limits;
│ │ │ │ +
799 for (
int i=0; i<dim; i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
802 if (currentKnotSpan[i]<
order_[i])
│ │ │ │ +
803 limits[i] -= (
order_[i] - currentKnotSpan[i]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
809 std::array<unsigned int, dim> offset;
│ │ │ │ +
810 for (
int i=0; i<dim; i++)
│ │ │ │ +
811 offset[i] = std::max((
int)(currentKnotSpan[i] -
order_[i]),0);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
814 std::array<std::vector<R>, dim> oneDValues;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
817 std::array<std::vector<R>, dim> lowOrderOneDValues;
│ │ │ │ +
│ │ │ │ +
819 std::array<DynamicMatrix<R>, dim> values;
│ │ │ │ +
│ │ │ │ +
821 for (
size_t i=0; i<dim; i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
825 for (
size_t j=0; j<oneDValues[i].size(); j++)
│ │ │ │ +
826 oneDValues[i][j] = values[i][
order_[i]][j];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
831 for (
size_t j=0; j<lowOrderOneDValues[i].size(); j++)
│ │ │ │ +
832 lowOrderOneDValues[i][j] = values[i][
order_[i]-1][j];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
838 std::array<std::vector<R>, dim> oneDDerivatives;
│ │ │ │ +
839 for (
size_t i=0; i<dim; i++)
│ │ │ │ +
│ │ │ │ +
841 oneDDerivatives[i].resize(limits[i]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
844 std::fill(oneDDerivatives[i].begin(), oneDDerivatives[i].end(),
R(0.0));
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
847 for (
size_t j=offset[i]; j<offset[i]+limits[i]; j++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
852 if (std::isnan(derivativeAddend1))
│ │ │ │ +
853 derivativeAddend1 = 0;
│ │ │ │ +
854 if (std::isnan(derivativeAddend2))
│ │ │ │ +
855 derivativeAddend2 = 0;
│ │ │ │ +
856 oneDDerivatives[i][j-offset[i]] =
order_[i] * ( derivativeAddend1 - derivativeAddend2 );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
863 std::array<std::vector<R>, dim> oneDValuesShort;
│ │ │ │ +
│ │ │ │ +
865 for (
int i=0; i<dim; i++)
│ │ │ │ +
│ │ │ │ +
867 oneDValuesShort[i].resize(limits[i]);
│ │ │ │ +
│ │ │ │ +
869 for (
size_t j=0; j<limits[i]; j++)
│ │ │ │ +
870 oneDValuesShort[i][j] = oneDValues[i][offset[i] + j];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
876 MultiDigitCounter ijkCounter(limits);
│ │ │ │ +
│ │ │ │ +
878 out.resize(ijkCounter.cycle());
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
881 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
│ │ │ │ +
882 for (
int j=0; j<dim; j++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
885 for (
int k=0; k<dim; k++)
│ │ │ │ +
886 out[i][0][j] *= (j==k) ? oneDDerivatives[k][ijkCounter[k]]
│ │ │ │ +
887 : oneDValuesShort[k][ijkCounter[k]];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
893 template <
size_type k>
│ │ │ │ +
│ │ │ │ +
894 void evaluate(
const typename std::array<int,k>& directions,
│ │ │ │ +
895 const FieldVector<typename GV::ctype,dim>& in,
│ │ │ │ +
896 std::vector<FieldVector<R,1> >& out,
│ │ │ │ +
897 const std::array<unsigned,dim>& currentKnotSpan)
const
│ │ │ │ +
│ │ │ │ +
899 if (k != 1 && k != 2)
│ │ │ │ +
900 DUNE_THROW(RangeError,
"Differentiation order greater than 2 is not supported!");
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
903 std::array<std::vector<R>, dim> oneDValues;
│ │ │ │ +
904 std::array<std::vector<R>, dim> oneDDerivatives;
│ │ │ │ +
905 std::array<std::vector<R>, dim> oneDSecondDerivatives;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
909 for (
size_t i=0; i<dim; i++)
│ │ │ │ +
910 evaluateAll(in[i], oneDValues[i],
true, oneDDerivatives[i],
false, oneDSecondDerivatives[i],
knotVectors_[i],
order_[i], currentKnotSpan[i]);
│ │ │ │ +
│ │ │ │ +
912 for (
size_t i=0; i<dim; i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
916 std::array<unsigned int, dim> offset;
│ │ │ │ +
917 for (
int i=0; i<dim; i++)
│ │ │ │ +
918 offset[i] = std::max((
int)(currentKnotSpan[i] -
order_[i]),0);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
921 std::array<unsigned int, dim> limits;
│ │ │ │ +
922 for (
int i=0; i<dim; i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
927 if (currentKnotSpan[i]<
order_[i])
│ │ │ │ +
928 limits[i] -= (
order_[i] - currentKnotSpan[i]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
935 std::array<std::vector<R>, dim> oneDValuesShort;
│ │ │ │ +
│ │ │ │ +
937 for (
int i=0; i<dim; i++)
│ │ │ │ +
│ │ │ │ +
939 oneDValuesShort[i].resize(limits[i]);
│ │ │ │ +
│ │ │ │ +
941 for (
size_t j=0; j<limits[i]; j++)
│ │ │ │ +
942 oneDValuesShort[i][j] = oneDValues[i][offset[i] + j];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
946 MultiDigitCounter ijkCounter(limits);
│ │ │ │ +
│ │ │ │ +
948 out.resize(ijkCounter.cycle());
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
953 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
956 for (
int l=0; l<dim; l++)
│ │ │ │ +
957 out[i][0] *= (directions[0]==l) ? oneDDerivatives[l][ijkCounter[l]]
│ │ │ │ +
958 : oneDValuesShort[l][ijkCounter[l]];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
965 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
968 for (
int j=0; j<dim; j++)
│ │ │ │ +
│ │ │ │ +
970 if (directions[0] != directions[1])
│ │ │ │ +
971 if (directions[0] == j || directions[1] == j)
│ │ │ │ +
972 out[i][0] *= oneDDerivatives[j][ijkCounter[j]];
│ │ │ │ +
│ │ │ │ +
974 out[i][0] *= oneDValuesShort[j][ijkCounter[j]];
│ │ │ │ +
│ │ │ │ +
976 if (directions[0] == j)
│ │ │ │ +
977 out[i][0] *= oneDSecondDerivatives[j][ijkCounter[j]];
│ │ │ │ +
│ │ │ │ +
979 out[i][0] *= oneDValuesShort[j][ijkCounter[j]];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
990 static std::array<unsigned int,dim>
getIJK(
typename GridView::IndexSet::IndexType idx, std::array<unsigned int,dim> elements)
│ │ │ │ +
│ │ │ │ +
992 std::array<unsigned,dim> result;
│ │ │ │ +
993 for (
int i=0; i<dim; i++)
│ │ │ │ +
│ │ │ │ +
995 result[i] = idx%elements[i];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1010 const std::vector<R>& knotVector,
│ │ │ │ +
│ │ │ │ +
1012 unsigned int currentKnotSpan)
│ │ │ │ +
│ │ │ │ +
1014 std::size_t outSize = order+1;
│ │ │ │ +
1015 if (currentKnotSpan<order)
│ │ │ │ +
1016 outSize -= (order - currentKnotSpan);
│ │ │ │ +
1017 if ( order > (knotVector.size() - currentKnotSpan - 2) )
│ │ │ │ +
1018 outSize -= order - (knotVector.size() - currentKnotSpan - 2);
│ │ │ │ +
1019 out.resize(outSize);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1022 DynamicMatrix<R> N(order+1, knotVector.size()-1);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1030 for (
size_t i=0; i<knotVector.size()-1; i++)
│ │ │ │ +
1031 N[0][i] = (i == currentKnotSpan);
│ │ │ │ +
│ │ │ │ +
1033 for (
size_t r=1; r<=order; r++)
│ │ │ │ +
1034 for (
size_t i=0; i<knotVector.size()-r-1; i++)
│ │ │ │ +
│ │ │ │ +
1036 R factor1 = ((knotVector[i+r] - knotVector[i]) > 1e-10)
│ │ │ │ +
1037 ? (in - knotVector[i]) / (knotVector[i+r] - knotVector[i])
│ │ │ │ +
│ │ │ │ +
1039 R factor2 = ((knotVector[i+r+1] - knotVector[i+1]) > 1e-10)
│ │ │ │ +
1040 ? (knotVector[i+r+1] - in) / (knotVector[i+r+1] - knotVector[i+1])
│ │ │ │ +
│ │ │ │ +
1042 N[r][i] = factor1 * N[r-1][i] + factor2 * N[r-1][i+1];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1049 for (
size_t i=0; i<out.size(); i++) {
│ │ │ │ +
1050 out[i] = N[order][std::max((
int)(currentKnotSpan - order),0) + i];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1067 DynamicMatrix<R>& out,
│ │ │ │ +
1068 const std::vector<R>& knotVector,
│ │ │ │ +
│ │ │ │ +
1070 unsigned int currentKnotSpan)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1073 DynamicMatrix<R>& N = out;
│ │ │ │ +
│ │ │ │ +
1075 N.resize(order+1, knotVector.size()-1);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1083 for (
size_t i=0; i<knotVector.size()-1; i++)
│ │ │ │ +
1084 N[0][i] = (i == currentKnotSpan);
│ │ │ │ +
│ │ │ │ +
1086 for (
size_t r=1; r<=order; r++)
│ │ │ │ +
1087 for (
size_t i=0; i<knotVector.size()-r-1; i++)
│ │ │ │ +
│ │ │ │ +
1089 R factor1 = ((knotVector[i+r] - knotVector[i]) > 1e-10)
│ │ │ │ +
1090 ? (in - knotVector[i]) / (knotVector[i+r] - knotVector[i])
│ │ │ │ +
│ │ │ │ +
1092 R factor2 = ((knotVector[i+r+1] - knotVector[i+1]) > 1e-10)
│ │ │ │ +
1093 ? (knotVector[i+r+1] - in) / (knotVector[i+r+1] - knotVector[i+1])
│ │ │ │ +
│ │ │ │ +
1095 N[r][i] = factor1 * N[r-1][i] + factor2 * N[r-1][i+1];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1110 std::vector<R>& out,
│ │ │ │ +
│ │ │ │ +
1112 bool evaluateHessian, std::vector<R>& outHess,
│ │ │ │ +
1113 const std::vector<R>& knotVector,
│ │ │ │ +
│ │ │ │ +
1115 unsigned int currentKnotSpan)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1120 if (currentKnotSpan<order)
│ │ │ │ +
1121 limit -= (order - currentKnotSpan);
│ │ │ │ +
1122 if ( order > (knotVector.size() - currentKnotSpan - 2) )
│ │ │ │ +
1123 limit -= order - (knotVector.size() - currentKnotSpan - 2);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1126 unsigned int offset;
│ │ │ │ +
1127 offset = std::max((
int)(currentKnotSpan - order),0);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1130 DynamicMatrix<R> values;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1134 out.resize(knotVector.size()-order-1);
│ │ │ │ +
1135 for (
size_t j=0; j<out.size(); j++)
│ │ │ │ +
1136 out[j] = values[order][j];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1139 std::vector<R> lowOrderOneDValues;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1143 lowOrderOneDValues.resize(knotVector.size()-(order-1)-1);
│ │ │ │ +
1144 for (
size_t j=0; j<lowOrderOneDValues.size(); j++)
│ │ │ │ +
1145 lowOrderOneDValues[j] = values[order-1][j];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1149 std::vector<R> lowOrderTwoDValues;
│ │ │ │ +
│ │ │ │ +
1151 if (order>1 && evaluateHessian)
│ │ │ │ +
│ │ │ │ +
1153 lowOrderTwoDValues.resize(knotVector.size()-(order-2)-1);
│ │ │ │ +
1154 for (
size_t j=0; j<lowOrderTwoDValues.size(); j++)
│ │ │ │ +
1155 lowOrderTwoDValues[j] = values[order-2][j];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1161 outJac.resize(limit);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1164 std::fill(outJac.begin(), outJac.end(),
R(0.0));
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1167 for (
size_t j=offset; j<offset+limit; j++)
│ │ │ │ +
│ │ │ │ +
1169 R derivativeAddend1 = lowOrderOneDValues[j] / (knotVector[j+order]-knotVector[j]);
│ │ │ │ +
1170 R derivativeAddend2 = lowOrderOneDValues[j+1] / (knotVector[j+order+1]-knotVector[j+1]);
│ │ │ │ +
│ │ │ │ +
1172 if (std::isnan(derivativeAddend1))
│ │ │ │ +
1173 derivativeAddend1 = 0;
│ │ │ │ +
1174 if (std::isnan(derivativeAddend2))
│ │ │ │ +
1175 derivativeAddend2 = 0;
│ │ │ │ +
1176 outJac[j-offset] = order * ( derivativeAddend1 - derivativeAddend2 );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1182 if (evaluateHessian)
│ │ │ │ +
│ │ │ │ +
1184 outHess.resize(limit);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1187 std::fill(outHess.begin(), outHess.end(),
R(0.0));
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1190 for (
size_t j=offset; j<offset+limit; j++)
│ │ │ │ +
│ │ │ │ +
1192 assert(j+2 < lowOrderTwoDValues.size());
│ │ │ │ +
1193 R derivativeAddend1 = lowOrderTwoDValues[j] / (knotVector[j+order]-knotVector[j]) / (knotVector[j+order-1]-knotVector[j]);
│ │ │ │ +
1194 R derivativeAddend2 = lowOrderTwoDValues[j+1] / (knotVector[j+order]-knotVector[j]) / (knotVector[j+order]-knotVector[j+1]);
│ │ │ │ +
1195 R derivativeAddend3 = lowOrderTwoDValues[j+1] / (knotVector[j+order+1]-knotVector[j+1]) / (knotVector[j+order]-knotVector[j+1]);
│ │ │ │ +
1196 R derivativeAddend4 = lowOrderTwoDValues[j+2] / (knotVector[j+order+1]-knotVector[j+1]) / (knotVector[j+1+order]-knotVector[j+2]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1199 if (std::isnan(derivativeAddend1))
│ │ │ │ +
1200 derivativeAddend1 = 0;
│ │ │ │ +
1201 if (std::isnan(derivativeAddend2))
│ │ │ │ +
1202 derivativeAddend2 = 0;
│ │ │ │ +
1203 if (std::isnan(derivativeAddend3))
│ │ │ │ +
1204 derivativeAddend3 = 0;
│ │ │ │ +
1205 if (std::isnan(derivativeAddend4))
│ │ │ │ +
1206 derivativeAddend4 = 0;
│ │ │ │ +
1207 outHess[j-offset] = order * (order-1) * ( derivativeAddend1 - derivativeAddend2 -derivativeAddend3 + derivativeAddend4 );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1228template<
typename GV>
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1232 static const int dim = GV::dimension;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1237 using Element =
typename GV::template Codim<0>::Entity;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1264 auto elementIndex =
preBasis_->gridView().indexSet().index(e);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1279namespace BasisFactory {
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1287inline auto bSpline(
const std::vector<double>& knotVector,
│ │ │ │ +
│ │ │ │ +
1289 bool makeOpen =
true)
│ │ │ │ +
│ │ │ │ +
1291 return [&knotVector, order, makeOpen](
const auto& gridView) {
│ │ │ │ +
1292 return BSplinePreBasis<std::decay_t<
decltype(gridView)>>(gridView, knotVector, order, makeOpen);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
1308template<
typename GV>
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
auto bSpline(const std::vector< double > &knotVector, unsigned int order, bool makeOpen=true)
Create a pre-basis factory that can create a B-spline pre-basis.
Definition bsplinebasis.hh:1287
│ │ │ │
Definition polynomial.hh:17
│ │ │ │ +
LocalFiniteElement in the sense of dune-localfunctions, for the B-spline basis on tensor-product grid...
Definition bsplinebasis.hh:367
│ │ │ │ +
BSplineLocalFiniteElement(const BSplineLocalFiniteElement &other)
Copy constructor.
Definition bsplinebasis.hh:388
│ │ │ │ +
const BSplinePreBasis< GV > & preBasis_
Definition bsplinebasis.hh:479
│ │ │ │ +
const BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R > > & localInterpolation() const
Hand out a LocalInterpolation object.
Definition bsplinebasis.hh:444
│ │ │ │ +
LocalFiniteElementTraits< BSplineLocalBasis< GV, R >, BSplineLocalCoefficients< dim >, BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R > > > Traits
Export various types related to this LocalFiniteElement.
Definition bsplinebasis.hh:377
│ │ │ │ +
std::array< unsigned, dim > currentKnotSpan_
Definition bsplinebasis.hh:486
│ │ │ │ +
BSplineLocalFiniteElement(const BSplinePreBasis< GV > &preBasis)
Constructor with a given B-spline basis.
Definition bsplinebasis.hh:381
│ │ │ │ +
const BSplineLocalCoefficients< dim > & localCoefficients() const
Hand out a LocalCoefficients object.
Definition bsplinebasis.hh:438
│ │ │ │ +
void bind(const std::array< unsigned, dim > &elementIdx)
Bind LocalFiniteElement to a specific knot span of the spline patch.
Definition bsplinebasis.hh:399
│ │ │ │ +
BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R > > localInterpolation_
Definition bsplinebasis.hh:483
│ │ │ │ +
GeometryType type() const
Return the reference element that the local finite element is defined on (here, a hypercube)
Definition bsplinebasis.hh:460
│ │ │ │ +
unsigned size() const
Number of shape functions in this finite element.
Definition bsplinebasis.hh:450
│ │ │ │ +
BSplineLocalCoefficients< dim > localCoefficients_
Definition bsplinebasis.hh:482
│ │ │ │ +
unsigned int size(int i) const
Number of degrees of freedom for one coordinate direction.
Definition bsplinebasis.hh:468
│ │ │ │ +
BSplineLocalBasis< GV, R > localBasis_
Definition bsplinebasis.hh:481
│ │ │ │ +
const BSplineLocalBasis< GV, R > & localBasis() const
Hand out a LocalBasis object.
Definition bsplinebasis.hh:432
│ │ │ │ +
Pre-basis for B-spline basis.
Definition bsplinebasis.hh:505
│ │ │ │ +
std::array< unsigned, dim > elements_
Number of grid elements in the different coordinate directions.
Definition bsplinebasis.hh:1221
│ │ │ │ +
GridView gridView_
Definition bsplinebasis.hh:1223
│ │ │ │ +
double R
Definition bsplinebasis.hh:574
│ │ │ │ +
static void evaluateFunctionFull(const typename GV::ctype &in, DynamicMatrix< R > &out, const std::vector< R > &knotVector, unsigned int order, unsigned int currentKnotSpan)
Evaluate all one-dimensional B-spline functions for a given coordinate direction.
Definition bsplinebasis.hh:1066
│ │ │ │ +
void evaluateFunction(const FieldVector< typename GV::ctype, dim > &in, std::vector< FieldVector< R, 1 > > &out, const std::array< unsigned, dim > ¤tKnotSpan) const
Evaluate all B-spline basis functions at a given point.
Definition bsplinebasis.hh:762
│ │ │ │ +
std::array< unsigned int, dim > order_
Order of the B-spline for each space dimension.
Definition bsplinebasis.hh:1215
│ │ │ │ +
static void evaluateAll(const typename GV::ctype &in, std::vector< R > &out, bool evaluateJacobian, std::vector< R > &outJac, bool evaluateHessian, std::vector< R > &outHess, const std::vector< R > &knotVector, unsigned int order, unsigned int currentKnotSpan)
Evaluate the second derivatives of all one-dimensional B-spline functions for a given coordinate dire...
Definition bsplinebasis.hh:1109
│ │ │ │ +
static void evaluateFunction(const typename GV::ctype &in, std::vector< R > &out, const std::vector< R > &knotVector, unsigned int order, unsigned int currentKnotSpan)
Evaluate all one-dimensional B-spline functions for a given coordinate direction.
Definition bsplinebasis.hh:1009
│ │ │ │ +
unsigned int size(size_t d) const
Number of shape functions in one direction.
Definition bsplinebasis.hh:753
│ │ │ │ +
GV GridView
The grid view that the FE space is defined on.
Definition bsplinebasis.hh:568
│ │ │ │ +
It indices(const Node &node, It it) const
Maps from subtree index set [0..size-1] to a globally unique multi index in global basis.
Definition bsplinebasis.hh:714
│ │ │ │ +
void evaluate(const typename std::array< int, k > &directions, const FieldVector< typename GV::ctype, dim > &in, std::vector< FieldVector< R, 1 > > &out, const std::array< unsigned, dim > ¤tKnotSpan) const
Evaluate Derivatives of all B-spline basis functions.
Definition bsplinebasis.hh:894
│ │ │ │ +
std::size_t size_type
Definition bsplinebasis.hh:569
│ │ │ │ +
void update(const GridView &gv)
Update the stored grid view, to be called if the grid has changed.
Definition bsplinebasis.hh:690
│ │ │ │ +
unsigned int dimension() const
Total number of B-spline basis functions.
Definition bsplinebasis.hh:744
│ │ │ │ +
void initializeIndices()
Initialize the global indices.
Definition bsplinebasis.hh:680
│ │ │ │ +
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition bsplinebasis.hh:684
│ │ │ │ +
static std::array< unsigned int, dim > getIJK(typename GridView::IndexSet::IndexType idx, std::array< unsigned int, dim > elements)
Compute integer element coordinates from the element index.
Definition bsplinebasis.hh:990
│ │ │ │ +
Node makeNode() const
Create tree node.
Definition bsplinebasis.hh:698
│ │ │ │ +
BSplinePreBasis(const GridView &gridView, const std::vector< double > &knotVector, unsigned int order, bool makeOpen=true)
Construct a B-spline basis for a given grid view and set of knot vectors.
Definition bsplinebasis.hh:594
│ │ │ │ +
void evaluateJacobian(const FieldVector< typename GV::ctype, dim > &in, std::vector< FieldMatrix< R, 1, dim > > &out, const std::array< unsigned, dim > ¤tKnotSpan) const
Evaluate Jacobian of all B-spline basis functions.
Definition bsplinebasis.hh:793
│ │ │ │ +
std::array< std::vector< double >, dim > knotVectors_
The knot vectors, one for each space dimension.
Definition bsplinebasis.hh:1218
│ │ │ │ +
size_type maxNodeSize() const
Get the maximal number of DOFs associated to node for any element.
Definition bsplinebasis.hh:704
│ │ │ │ +
BSplinePreBasis(const GridView &gridView, const FieldVector< double, dim > &lowerLeft, const FieldVector< double, dim > &upperRight, const std::array< unsigned int, dim > &elements, unsigned int order, bool makeOpen=true)
Construct a B-spline basis for a given grid view with uniform knot vectors.
Definition bsplinebasis.hh:646
│ │ │ │ +
LocalBasis class in the sense of dune-localfunctions, presenting the restriction of a B-spline patch ...
Definition bsplinebasis.hh:51
│ │ │ │ +
LocalBasisTraits< D, dim, FieldVector< D, dim >, R, 1, FieldVector< R, 1 >, FieldMatrix< R, 1, dim > > Traits
export type traits for function signature
Definition bsplinebasis.hh:60
│ │ │ │ +
unsigned int order() const
Polynomial order of the shape functions.
Definition bsplinebasis.hh:145
│ │ │ │ +
std::size_t size() const
Return the number of basis functions on the current knot span.
Definition bsplinebasis.hh:152
│ │ │ │ +
void evaluate(const typename std::array< int, k > &directions, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions and derivatives of any order.
Definition bsplinebasis.hh:102
│ │ │ │ +
void evaluateFunction(const FieldVector< D, dim > &in, std::vector< FieldVector< R, 1 > > &out) const
Evaluate all shape functions.
Definition bsplinebasis.hh:75
│ │ │ │ +
void evaluateJacobian(const FieldVector< D, dim > &in, std::vector< FieldMatrix< D, 1, dim > > &out) const
Evaluate Jacobian of all shape functions.
Definition bsplinebasis.hh:87
│ │ │ │ +
BSplineLocalBasis(const BSplinePreBasis< GV > &preBasis, const BSplineLocalFiniteElement< GV, R > &lFE)
Constructor with a given B-spline patch.
Definition bsplinebasis.hh:66
│ │ │ │ +
Attaches a shape function to an entity.
Definition bsplinebasis.hh:183
│ │ │ │ +
const LocalKey & localKey(std::size_t i) const
get i'th index
Definition bsplinebasis.hh:326
│ │ │ │ +
void init(const std::array< unsigned, dim > &sizes)
Definition bsplinebasis.hh:261
│ │ │ │ +
std::size_t size() const
number of coefficients
Definition bsplinebasis.hh:320
│ │ │ │ +
Local interpolation in the sense of dune-localfunctions, for the B-spline basis on tensor-product gri...
Definition bsplinebasis.hh:345
│ │ │ │ +
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition bsplinebasis.hh:349
│ │ │ │ +
Definition bsplinebasis.hh:1231
│ │ │ │ +
const FiniteElement & finiteElement() const
Return the LocalFiniteElement for the element we are bound to.
Definition bsplinebasis.hh:1255
│ │ │ │ +
typename GV::template Codim< 0 >::Entity Element
Definition bsplinebasis.hh:1237
│ │ │ │ +
const BSplinePreBasis< GV > * preBasis_
Definition bsplinebasis.hh:1271
│ │ │ │ +
Element element_
Definition bsplinebasis.hh:1274
│ │ │ │ +
void bind(const Element &e)
Bind to element.
Definition bsplinebasis.hh:1261
│ │ │ │ +
BSplineNode(const BSplinePreBasis< GV > *preBasis)
Definition bsplinebasis.hh:1240
│ │ │ │ +
const Element & element() const
Return current element, throw if unbound.
Definition bsplinebasis.hh:1246
│ │ │ │ +
FiniteElement finiteElement_
Definition bsplinebasis.hh:1273
│ │ │ │ +
std::size_t size_type
Definition bsplinebasis.hh:1236
│ │ │ │ +
Global basis for given pre-basis.
Definition defaultglobalbasis.hh:50
│ │ │ │ +
A generic MixIn class for PreBasis.
Definition leafprebasismixin.hh:36
│ │ │ │ +
size_type size() const
Get the total dimension of the space spanned by this basis.
Definition leafprebasismixin.hh:60
│ │ │ │ +
size_type size() const
Definition nodes.hh:147
│ │ │ │ +
void setSize(const size_type size)
Definition nodes.hh:169
│ │ │ │ +
│ │ │ │
