│ │ │ │
│ │ │ │ -
3#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH
│ │ │ │ -
4#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BSPLINEBASIS_HH
│ │ │ │ +
3#ifndef DUNE_FUNCTIONS_GRIDFUNCTIONS_GRIDVIEWENTITYSET_HH
│ │ │ │ +
4#define DUNE_FUNCTIONS_GRIDFUNCTIONS_GRIDVIEWENTITYSET_HH
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
14#include <dune/common/dynmatrix.hh>
│ │ │ │ -
│ │ │ │ -
16#include <dune/localfunctions/common/localbasis.hh>
│ │ │ │ -
17#include <dune/common/diagonalmatrix.hh>
│ │ │ │ -
18#include <dune/localfunctions/common/localkey.hh>
│ │ │ │ -
19#include <dune/localfunctions/common/localfiniteelementtraits.hh>
│ │ │ │ -
20#include <dune/geometry/type.hh>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
29template<
typename GV,
typename R>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
44template<
class GV,
class R>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
49 typedef typename GV::ctype D;
│ │ │ │ -
50 enum {dim = GV::dimension};
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
54 typedef LocalBasisTraits<D,dim,FieldVector<D,dim>,R,1,FieldVector<R,1>,
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
21template<
class GV,
int cd>
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
32 typedef typename GridView::template Codim<codim>::Entity
Element;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
54 return gv_.contains(e);
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
63 : preBasis_(preBasis),
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
60 return gv_.size(
codim);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
66 return gv_.template begin<codim>();
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
71 std::vector<FieldVector<R,1> >& out)
const
│ │ │ │ -
│ │ │ │ -
73 FieldVector<D,dim> globalIn = offset_;
│ │ │ │ -
74 scaling_.umv(in,globalIn);
│ │ │ │ -
│ │ │ │ -
76 preBasis_.evaluateFunction(globalIn, out, lFE_.currentKnotSpan_);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
83 std::vector<FieldMatrix<D,1,dim> >& out)
const
│ │ │ │ -
│ │ │ │ -
85 FieldVector<D,dim> globalIn = offset_;
│ │ │ │ -
86 scaling_.umv(in,globalIn);
│ │ │ │ -
│ │ │ │ -
88 preBasis_.evaluateJacobian(globalIn, out, lFE_.currentKnotSpan_);
│ │ │ │ -
│ │ │ │ -
90 for (
size_t i=0; i<out.size(); i++)
│ │ │ │ -
91 for (
int j=0; j<dim; j++)
│ │ │ │ -
92 out[i][0][j] *= scaling_[j][j];
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
97 inline void evaluate (
const typename std::array<int,k>& directions,
│ │ │ │ -
98 const typename Traits::DomainType& in,
│ │ │ │ -
99 std::vector<typename Traits::RangeType>& out)
const
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
108 FieldVector<D,dim> globalIn = offset_;
│ │ │ │ -
109 scaling_.umv(in,globalIn);
│ │ │ │ -
│ │ │ │ -
111 preBasis_.evaluate(directions, globalIn, out, lFE_.currentKnotSpan_);
│ │ │ │ -
│ │ │ │ -
113 for (
size_t i=0; i<out.size(); i++)
│ │ │ │ -
114 out[i][0] *= scaling_[directions[0]][directions[0]];
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
119 FieldVector<D,dim> globalIn = offset_;
│ │ │ │ -
120 scaling_.umv(in,globalIn);
│ │ │ │ -
│ │ │ │ -
122 preBasis_.evaluate(directions, globalIn, out, lFE_.currentKnotSpan_);
│ │ │ │ -
│ │ │ │ -
124 for (
size_t i=0; i<out.size(); i++)
│ │ │ │ -
125 out[i][0] *= scaling_[directions[0]][directions[0]]*scaling_[directions[1]][directions[1]];
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
129 DUNE_THROW(NotImplemented,
"B-Spline derivatives of order " << k <<
" not implemented yet!");
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
142 return *std::max_element(preBasis_.order_.begin(), preBasis_.order_.end());
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
159 FieldVector<D,dim> offset_;
│ │ │ │ -
160 DiagonalMatrix<D,dim> scaling_;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
180 std::array<unsigned int,dim> multiindex (
unsigned int i)
const
│ │ │ │ -
│ │ │ │ -
182 std::array<unsigned int,dim> alpha;
│ │ │ │ -
183 for (
int j=0; j<dim; j++)
│ │ │ │ -
│ │ │ │ -
185 alpha[j] = i % sizes_[j];
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
192 void setup1d(std::vector<unsigned int>& subEntity)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
203 unsigned lastIndex=0;
│ │ │ │ -
204 subEntity[lastIndex++] = 0;
│ │ │ │ -
205 for (
unsigned i = 0; i < sizes_[0] - 2; ++i)
│ │ │ │ -
206 subEntity[lastIndex++] = 0;
│ │ │ │ -
│ │ │ │ -
208 subEntity[lastIndex++] = 1;
│ │ │ │ -
│ │ │ │ -
210 assert(
size()==lastIndex);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
213 void setup2d(std::vector<unsigned int>& subEntity)
│ │ │ │ -
│ │ │ │ -
215 unsigned lastIndex=0;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
229 subEntity[lastIndex++] = 0;
│ │ │ │ -
230 for (
unsigned i = 0; i < sizes_[0]-2; ++i)
│ │ │ │ -
231 subEntity[lastIndex++] = 2;
│ │ │ │ -
│ │ │ │ -
233 subEntity[lastIndex++] = 1;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
236 for (
unsigned e = 0; e < sizes_[1]-2; ++e)
│ │ │ │ -
│ │ │ │ -
238 subEntity[lastIndex++] = 0;
│ │ │ │ -
239 for (
unsigned i = 0; i < sizes_[0]-2; ++i)
│ │ │ │ -
240 subEntity[lastIndex++] = 0;
│ │ │ │ -
241 subEntity[lastIndex++] = 1;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
245 subEntity[lastIndex++] = 2;
│ │ │ │ -
246 for (
unsigned i = 0; i < sizes_[0]-2; ++i)
│ │ │ │ -
247 subEntity[lastIndex++] = 3;
│ │ │ │ -
│ │ │ │ -
249 subEntity[lastIndex++] = 3;
│ │ │ │ -
│ │ │ │ -
251 assert(
size()==lastIndex);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
256 void init(
const std::array<unsigned,dim>& sizes)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
263 std::vector<unsigned int> codim(li_.size());
│ │ │ │ -
│ │ │ │ -
265 for (std::size_t i=0; i<codim.size(); i++)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
270 std::array<unsigned int,dim> mIdx = multiindex(i);
│ │ │ │ -
271 for (
int j=0; j<dim; j++)
│ │ │ │ -
272 if (mIdx[j]==0 or mIdx[j]==sizes[j]-1)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
281 std::vector<unsigned int> index(
size());
│ │ │ │ -
│ │ │ │ -
283 for (std::size_t i=0; i<index.size(); i++)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
287 std::array<unsigned int,dim> mIdx = multiindex(i);
│ │ │ │ -
│ │ │ │ -
289 for (
int j=dim-1; j>=0; j--)
│ │ │ │ -
290 if (mIdx[j]>0 and mIdx[j]<sizes[j]-1)
│ │ │ │ -
291 index[i] = (sizes[j]-1)*index[i] + (mIdx[j]-1);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
295 std::vector<unsigned int> subEntity(li_.size());
│ │ │ │ -
│ │ │ │ -
297 if (subEntity.size() > 0)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
303 }
else if (dim==2 and sizes_[0]>1 and sizes_[1]>1) {
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
310 for (
size_t i=0; i<li_.size(); i++)
│ │ │ │ -
311 li_[i] = LocalKey(subEntity[i], codim[i], index[i]);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
317 return std::accumulate(sizes_.begin(), sizes_.end(), 1, std::multiplies<unsigned int>());
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
329 std::array<unsigned, dim> sizes_;
│ │ │ │ -
│ │ │ │ -
331 std::vector<LocalKey> li_;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
338template<
int dim,
class LB>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
343 template<
typename F,
typename C>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
346 DUNE_THROW(NotImplemented,
"BSplineLocalInterpolation::interpolate");
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
360template<
class GV,
class R>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
363 typedef typename GV::ctype D;
│ │ │ │ -
364 enum {dim = GV::dimension};
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
370 typedef LocalFiniteElementTraits<BSplineLocalBasis<GV,R>,
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
394 void bind(
const std::array<unsigned,dim>& elementIdx)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
397 for (
size_t i=0; i<elementIdx.size(); i++)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
405 for (
size_t j=0; j<elementIdx[i]; j++)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
420 std::array<unsigned int, dim> sizes;
│ │ │ │ -
421 for (
size_t i=0; i<dim; i++)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
448 for (
int i=0; i<dim; i++)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
457 return GeometryTypes::cube(dim);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
466 unsigned int r = order[i]+1;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
500 static const int dim = GV::dimension;
│ │ │ │ -
│ │ │ │ -
503 class MultiDigitCounter
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
510 MultiDigitCounter(
const std::array<unsigned int,dim>& limits)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
513 std::fill(counter_.begin(), counter_.end(), 0);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
517 MultiDigitCounter& operator++()
│ │ │ │ -
│ │ │ │ -
519 for (
int i=0; i<dim; i++)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
524 if (counter_[i] < limits_[i])
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
533 const unsigned int& operator[](
int i)
const
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
539 unsigned int cycle()
const
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
542 for (
int i=0; i<dim; i++)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
550 const std::array<unsigned int,dim> limits_;
│ │ │ │ -
│ │ │ │ -
553 std::array<unsigned int,dim> counter_;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
591 const std::vector<double>& knotVector,
│ │ │ │ -
│ │ │ │ -
593 bool makeOpen =
true)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
603 for (
int i=0; i<dim; i++)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
608 for (
unsigned int j=0; j<order; j++)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
614 for (
unsigned int j=0; j<order; j++)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
643 const FieldVector<double,dim>& lowerLeft,
│ │ │ │ -
644 const FieldVector<double,dim>& upperRight,
│ │ │ │ -
645 const std::array<unsigned int,dim>& elements,
│ │ │ │ -
│ │ │ │ -
647 bool makeOpen =
true)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
655 for (
int i=0; i<dim; i++)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
660 for (
unsigned int j=0; j<order; j++)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
664 for (
size_t j=0; j<elements[i]+1; j++)
│ │ │ │ -
665 knotVectors_[i].push_back(lowerLeft[i] + j*(upperRight[i]-lowerLeft[i]) / elements[i]);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
668 for (
unsigned int j=0; j<order; j++)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
711 template<
class ST,
int i>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
714 assert(prefix.size() == 0 || prefix.size() == 1);
│ │ │ │ -
715 return (prefix.size() == 0) ?
size() : 0;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
728 for (
int i=0; i<dim; i++)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
734 template<
typename It>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
739 std::array<unsigned int, dim> localSizes;
│ │ │ │ -
740 for (
int i=0; i<dim; i++)
│ │ │ │ -
│ │ │ │ -
742 for (
size_type i = 0, end = node.
size() ; i < end ; ++i, ++it)
│ │ │ │ -
│ │ │ │ -
744 std::array<unsigned int,dim> localIJK =
getIJK(i, localSizes);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
747 const auto order =
order_;
│ │ │ │ -
│ │ │ │ -
749 std::array<unsigned int,dim> globalIJK;
│ │ │ │ -
750 for (
int i=0; i<dim; i++)
│ │ │ │ -
751 globalIJK[i] = std::max((
int)currentKnotSpan[i] - (
int)order[i], 0) + localIJK[i];
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
756 for (
int i=dim-2; i>=0; i--)
│ │ │ │ -
757 globalIdx = globalIdx *
size(i) + globalIJK[i];
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
767 unsigned int result = 1;
│ │ │ │ -
768 for (
size_t i=0; i<dim; i++)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
774 unsigned int size (
size_t d)
const
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
782 std::vector<FieldVector<R,1> >& out,
│ │ │ │ -
783 const std::array<unsigned,dim>& currentKnotSpan)
const
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
786 std::array<std::vector<R>, dim> oneDValues;
│ │ │ │ -
│ │ │ │ -
788 for (
size_t i=0; i<dim; i++)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
791 std::array<unsigned int, dim> limits;
│ │ │ │ -
792 for (
int i=0; i<dim; i++)
│ │ │ │ -
793 limits[i] = oneDValues[i].
size();
│ │ │ │ -
│ │ │ │ -
795 MultiDigitCounter ijkCounter(limits);
│ │ │ │ -
│ │ │ │ -
797 out.resize(ijkCounter.cycle());
│ │ │ │ -
│ │ │ │ -
799 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
802 for (
size_t j=0; j<dim; j++)
│ │ │ │ -
803 out[i] *= oneDValues[j][ijkCounter[j]];
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
813 std::vector<FieldMatrix<R,1,dim> >& out,
│ │ │ │ -
814 const std::array<unsigned,dim>& currentKnotSpan)
const
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
817 std::array<unsigned int, dim> limits;
│ │ │ │ -
818 for (
int i=0; i<dim; i++)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
821 if (currentKnotSpan[i]<
order_[i])
│ │ │ │ -
822 limits[i] -= (
order_[i] - currentKnotSpan[i]);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
828 std::array<unsigned int, dim> offset;
│ │ │ │ -
829 for (
int i=0; i<dim; i++)
│ │ │ │ -
830 offset[i] = std::max((
int)(currentKnotSpan[i] -
order_[i]),0);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
833 std::array<std::vector<R>, dim> oneDValues;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
836 std::array<std::vector<R>, dim> lowOrderOneDValues;
│ │ │ │ -
│ │ │ │ -
838 std::array<DynamicMatrix<R>, dim> values;
│ │ │ │ -
│ │ │ │ -
840 for (
size_t i=0; i<dim; i++)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
844 for (
size_t j=0; j<oneDValues[i].size(); j++)
│ │ │ │ -
845 oneDValues[i][j] = values[i][
order_[i]][j];
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
850 for (
size_t j=0; j<lowOrderOneDValues[i].size(); j++)
│ │ │ │ -
851 lowOrderOneDValues[i][j] = values[i][
order_[i]-1][j];
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
857 std::array<std::vector<R>, dim> oneDDerivatives;
│ │ │ │ -
858 for (
size_t i=0; i<dim; i++)
│ │ │ │ -
│ │ │ │ -
860 oneDDerivatives[i].resize(limits[i]);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
863 std::fill(oneDDerivatives[i].begin(), oneDDerivatives[i].end(),
R(0.0));
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
866 for (
size_t j=offset[i]; j<offset[i]+limits[i]; j++)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
871 if (std::isnan(derivativeAddend1))
│ │ │ │ -
872 derivativeAddend1 = 0;
│ │ │ │ -
873 if (std::isnan(derivativeAddend2))
│ │ │ │ -
874 derivativeAddend2 = 0;
│ │ │ │ -
875 oneDDerivatives[i][j-offset[i]] =
order_[i] * ( derivativeAddend1 - derivativeAddend2 );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
882 std::array<std::vector<R>, dim> oneDValuesShort;
│ │ │ │ -
│ │ │ │ -
884 for (
int i=0; i<dim; i++)
│ │ │ │ -
│ │ │ │ -
886 oneDValuesShort[i].resize(limits[i]);
│ │ │ │ -
│ │ │ │ -
888 for (
size_t j=0; j<limits[i]; j++)
│ │ │ │ -
889 oneDValuesShort[i][j] = oneDValues[i][offset[i] + j];
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
895 MultiDigitCounter ijkCounter(limits);
│ │ │ │ -
│ │ │ │ -
897 out.resize(ijkCounter.cycle());
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
900 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
│ │ │ │ -
901 for (
int j=0; j<dim; j++)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
904 for (
int k=0; k<dim; k++)
│ │ │ │ -
905 out[i][0][j] *= (j==k) ? oneDDerivatives[k][ijkCounter[k]]
│ │ │ │ -
906 : oneDValuesShort[k][ijkCounter[k]];
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
912 template <
size_type k>
│ │ │ │ -
│ │ │ │ -
913 void evaluate(
const typename std::array<int,k>& directions,
│ │ │ │ -
914 const FieldVector<typename GV::ctype,dim>& in,
│ │ │ │ -
915 std::vector<FieldVector<R,1> >& out,
│ │ │ │ -
916 const std::array<unsigned,dim>& currentKnotSpan)
const
│ │ │ │ -
│ │ │ │ -
918 if (k != 1 && k != 2)
│ │ │ │ -
919 DUNE_THROW(RangeError,
"Differentiation order greater than 2 is not supported!");
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
922 std::array<std::vector<R>, dim> oneDValues;
│ │ │ │ -
923 std::array<std::vector<R>, dim> oneDDerivatives;
│ │ │ │ -
924 std::array<std::vector<R>, dim> oneDSecondDerivatives;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
928 for (
size_t i=0; i<dim; i++)
│ │ │ │ -
929 evaluateAll(in[i], oneDValues[i],
true, oneDDerivatives[i],
false, oneDSecondDerivatives[i],
knotVectors_[i],
order_[i], currentKnotSpan[i]);
│ │ │ │ -
│ │ │ │ -
931 for (
size_t i=0; i<dim; i++)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
935 std::array<unsigned int, dim> offset;
│ │ │ │ -
936 for (
int i=0; i<dim; i++)
│ │ │ │ -
937 offset[i] = std::max((
int)(currentKnotSpan[i] -
order_[i]),0);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
940 std::array<unsigned int, dim> limits;
│ │ │ │ -
941 for (
int i=0; i<dim; i++)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
946 if (currentKnotSpan[i]<
order_[i])
│ │ │ │ -
947 limits[i] -= (
order_[i] - currentKnotSpan[i]);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
954 std::array<std::vector<R>, dim> oneDValuesShort;
│ │ │ │ -
│ │ │ │ -
956 for (
int i=0; i<dim; i++)
│ │ │ │ -
│ │ │ │ -
958 oneDValuesShort[i].resize(limits[i]);
│ │ │ │ -
│ │ │ │ -
960 for (
size_t j=0; j<limits[i]; j++)
│ │ │ │ -
961 oneDValuesShort[i][j] = oneDValues[i][offset[i] + j];
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
965 MultiDigitCounter ijkCounter(limits);
│ │ │ │ -
│ │ │ │ -
967 out.resize(ijkCounter.cycle());
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
972 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
975 for (
int l=0; l<dim; l++)
│ │ │ │ -
976 out[i][0] *= (directions[0]==l) ? oneDDerivatives[l][ijkCounter[l]]
│ │ │ │ -
977 : oneDValuesShort[l][ijkCounter[l]];
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
984 for (
size_t i=0; i<out.size(); i++, ++ijkCounter)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
987 for (
int j=0; j<dim; j++)
│ │ │ │ -
│ │ │ │ -
989 if (directions[0] != directions[1])
│ │ │ │ -
990 if (directions[0] == j || directions[1] == j)
│ │ │ │ -
991 out[i][0] *= oneDDerivatives[j][ijkCounter[j]];
│ │ │ │ -
│ │ │ │ -
993 out[i][0] *= oneDValuesShort[j][ijkCounter[j]];
│ │ │ │ -
│ │ │ │ -
995 if (directions[0] == j)
│ │ │ │ -
996 out[i][0] *= oneDSecondDerivatives[j][ijkCounter[j]];
│ │ │ │ -
│ │ │ │ -
998 out[i][0] *= oneDValuesShort[j][ijkCounter[j]];
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1009 static std::array<unsigned int,dim>
getIJK(
typename GridView::IndexSet::IndexType idx, std::array<unsigned int,dim> elements)
│ │ │ │ -
│ │ │ │ -
1011 std::array<unsigned,dim> result;
│ │ │ │ -
1012 for (
int i=0; i<dim; i++)
│ │ │ │ -
│ │ │ │ -
1014 result[i] = idx%elements[i];
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1029 const std::vector<R>& knotVector,
│ │ │ │ -
│ │ │ │ -
1031 unsigned int currentKnotSpan)
│ │ │ │ -
│ │ │ │ -
1033 std::size_t outSize = order+1;
│ │ │ │ -
1034 if (currentKnotSpan<order)
│ │ │ │ -
1035 outSize -= (order - currentKnotSpan);
│ │ │ │ -
1036 if ( order > (knotVector.size() - currentKnotSpan - 2) )
│ │ │ │ -
1037 outSize -= order - (knotVector.size() - currentKnotSpan - 2);
│ │ │ │ -
1038 out.resize(outSize);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1041 DynamicMatrix<R> N(order+1, knotVector.size()-1);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1049 for (
size_t i=0; i<knotVector.size()-1; i++)
│ │ │ │ -
1050 N[0][i] = (i == currentKnotSpan);
│ │ │ │ -
│ │ │ │ -
1052 for (
size_t r=1; r<=order; r++)
│ │ │ │ -
1053 for (
size_t i=0; i<knotVector.size()-r-1; i++)
│ │ │ │ -
│ │ │ │ -
1055 R factor1 = ((knotVector[i+r] - knotVector[i]) > 1e-10)
│ │ │ │ -
1056 ? (in - knotVector[i]) / (knotVector[i+r] - knotVector[i])
│ │ │ │ -
│ │ │ │ -
1058 R factor2 = ((knotVector[i+r+1] - knotVector[i+1]) > 1e-10)
│ │ │ │ -
1059 ? (knotVector[i+r+1] - in) / (knotVector[i+r+1] - knotVector[i+1])
│ │ │ │ -
│ │ │ │ -
1061 N[r][i] = factor1 * N[r-1][i] + factor2 * N[r-1][i+1];
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1068 for (
size_t i=0; i<out.size(); i++) {
│ │ │ │ -
1069 out[i] = N[order][std::max((
int)(currentKnotSpan - order),0) + i];
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1086 DynamicMatrix<R>& out,
│ │ │ │ -
1087 const std::vector<R>& knotVector,
│ │ │ │ -
│ │ │ │ -
1089 unsigned int currentKnotSpan)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1092 DynamicMatrix<R>& N = out;
│ │ │ │ -
│ │ │ │ -
1094 N.resize(order+1, knotVector.size()-1);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1102 for (
size_t i=0; i<knotVector.size()-1; i++)
│ │ │ │ -
1103 N[0][i] = (i == currentKnotSpan);
│ │ │ │ -
│ │ │ │ -
1105 for (
size_t r=1; r<=order; r++)
│ │ │ │ -
1106 for (
size_t i=0; i<knotVector.size()-r-1; i++)
│ │ │ │ -
│ │ │ │ -
1108 R factor1 = ((knotVector[i+r] - knotVector[i]) > 1e-10)
│ │ │ │ -
1109 ? (in - knotVector[i]) / (knotVector[i+r] - knotVector[i])
│ │ │ │ -
│ │ │ │ -
1111 R factor2 = ((knotVector[i+r+1] - knotVector[i+1]) > 1e-10)
│ │ │ │ -
1112 ? (knotVector[i+r+1] - in) / (knotVector[i+r+1] - knotVector[i+1])
│ │ │ │ -
│ │ │ │ -
1114 N[r][i] = factor1 * N[r-1][i] + factor2 * N[r-1][i+1];
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1129 std::vector<R>& out,
│ │ │ │ -
│ │ │ │ -
1131 bool evaluateHessian, std::vector<R>& outHess,
│ │ │ │ -
1132 const std::vector<R>& knotVector,
│ │ │ │ -
│ │ │ │ -
1134 unsigned int currentKnotSpan)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1139 if (currentKnotSpan<order)
│ │ │ │ -
1140 limit -= (order - currentKnotSpan);
│ │ │ │ -
1141 if ( order > (knotVector.size() - currentKnotSpan - 2) )
│ │ │ │ -
1142 limit -= order - (knotVector.size() - currentKnotSpan - 2);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1145 unsigned int offset;
│ │ │ │ -
1146 offset = std::max((
int)(currentKnotSpan - order),0);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1149 DynamicMatrix<R> values;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1153 out.resize(knotVector.size()-order-1);
│ │ │ │ -
1154 for (
size_t j=0; j<out.size(); j++)
│ │ │ │ -
1155 out[j] = values[order][j];
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1158 std::vector<R> lowOrderOneDValues;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1162 lowOrderOneDValues.resize(knotVector.size()-(order-1)-1);
│ │ │ │ -
1163 for (
size_t j=0; j<lowOrderOneDValues.size(); j++)
│ │ │ │ -
1164 lowOrderOneDValues[j] = values[order-1][j];
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1168 std::vector<R> lowOrderTwoDValues;
│ │ │ │ -
│ │ │ │ -
1170 if (order>1 && evaluateHessian)
│ │ │ │ -
│ │ │ │ -
1172 lowOrderTwoDValues.resize(knotVector.size()-(order-2)-1);
│ │ │ │ -
1173 for (
size_t j=0; j<lowOrderTwoDValues.size(); j++)
│ │ │ │ -
1174 lowOrderTwoDValues[j] = values[order-2][j];
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1180 outJac.resize(limit);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1183 std::fill(outJac.begin(), outJac.end(),
R(0.0));
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1186 for (
size_t j=offset; j<offset+limit; j++)
│ │ │ │ -
│ │ │ │ -
1188 R derivativeAddend1 = lowOrderOneDValues[j] / (knotVector[j+order]-knotVector[j]);
│ │ │ │ -
1189 R derivativeAddend2 = lowOrderOneDValues[j+1] / (knotVector[j+order+1]-knotVector[j+1]);
│ │ │ │ -
│ │ │ │ -
1191 if (std::isnan(derivativeAddend1))
│ │ │ │ -
1192 derivativeAddend1 = 0;
│ │ │ │ -
1193 if (std::isnan(derivativeAddend2))
│ │ │ │ -
1194 derivativeAddend2 = 0;
│ │ │ │ -
1195 outJac[j-offset] = order * ( derivativeAddend1 - derivativeAddend2 );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1201 if (evaluateHessian)
│ │ │ │ -
│ │ │ │ -
1203 outHess.resize(limit);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1206 std::fill(outHess.begin(), outHess.end(),
R(0.0));
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1209 for (
size_t j=offset; j<offset+limit; j++)
│ │ │ │ -
│ │ │ │ -
1211 assert(j+2 < lowOrderTwoDValues.size());
│ │ │ │ -
1212 R derivativeAddend1 = lowOrderTwoDValues[j] / (knotVector[j+order]-knotVector[j]) / (knotVector[j+order-1]-knotVector[j]);
│ │ │ │ -
1213 R derivativeAddend2 = lowOrderTwoDValues[j+1] / (knotVector[j+order]-knotVector[j]) / (knotVector[j+order]-knotVector[j+1]);
│ │ │ │ -
1214 R derivativeAddend3 = lowOrderTwoDValues[j+1] / (knotVector[j+order+1]-knotVector[j+1]) / (knotVector[j+order]-knotVector[j+1]);
│ │ │ │ -
1215 R derivativeAddend4 = lowOrderTwoDValues[j+2] / (knotVector[j+order+1]-knotVector[j+1]) / (knotVector[j+1+order]-knotVector[j+2]);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1218 if (std::isnan(derivativeAddend1))
│ │ │ │ -
1219 derivativeAddend1 = 0;
│ │ │ │ -
1220 if (std::isnan(derivativeAddend2))
│ │ │ │ -
1221 derivativeAddend2 = 0;
│ │ │ │ -
1222 if (std::isnan(derivativeAddend3))
│ │ │ │ -
1223 derivativeAddend3 = 0;
│ │ │ │ -
1224 if (std::isnan(derivativeAddend4))
│ │ │ │ -
1225 derivativeAddend4 = 0;
│ │ │ │ -
1226 outHess[j-offset] = order * (order-1) * ( derivativeAddend1 - derivativeAddend2 -derivativeAddend3 + derivativeAddend4 );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1247template<
typename GV>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1251 static const int dim = GV::dimension;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1256 using Element =
typename GV::template Codim<0>::Entity;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1283 auto elementIndex =
preBasis_->gridView().indexSet().index(e);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1298namespace BasisFactory {
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1306inline auto bSpline(
const std::vector<double>& knotVector,
│ │ │ │ -
│ │ │ │ -
1308 bool makeOpen =
true)
│ │ │ │ -
│ │ │ │ -
1310 return [&knotVector, order, makeOpen](
const auto& gridView) {
│ │ │ │ -
1311 return BSplinePreBasis<std::decay_t<
decltype(gridView)>>(gridView, knotVector, order, makeOpen);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
1327template<
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:1306
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
72 return gv_.template end<codim>();
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
Definition polynomial.hh:10
│ │ │ │ -
LocalFiniteElement in the sense of dune-localfunctions, for the B-spline basis on tensor-product grid...
Definition bsplinebasis.hh:362
│ │ │ │ -
BSplineLocalFiniteElement(const BSplineLocalFiniteElement &other)
Copy constructor.
Definition bsplinebasis.hh:383
│ │ │ │ -
const BSplinePreBasis< GV > & preBasis_
Definition bsplinebasis.hh:474
│ │ │ │ -
const BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R > > & localInterpolation() const
Hand out a LocalInterpolation object.
Definition bsplinebasis.hh:439
│ │ │ │ -
LocalFiniteElementTraits< BSplineLocalBasis< GV, R >, BSplineLocalCoefficients< dim >, BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R > > > Traits
Export various types related to this LocalFiniteElement.
Definition bsplinebasis.hh:372
│ │ │ │ -
std::array< unsigned, dim > currentKnotSpan_
Definition bsplinebasis.hh:481
│ │ │ │ -
BSplineLocalFiniteElement(const BSplinePreBasis< GV > &preBasis)
Constructor with a given B-spline basis.
Definition bsplinebasis.hh:376
│ │ │ │ -
const BSplineLocalCoefficients< dim > & localCoefficients() const
Hand out a LocalCoefficients object.
Definition bsplinebasis.hh:433
│ │ │ │ -
void bind(const std::array< unsigned, dim > &elementIdx)
Bind LocalFiniteElement to a specific knot span of the spline patch.
Definition bsplinebasis.hh:394
│ │ │ │ -
BSplineLocalInterpolation< dim, BSplineLocalBasis< GV, R > > localInterpolation_
Definition bsplinebasis.hh:478
│ │ │ │ -
GeometryType type() const
Return the reference element that the local finite element is defined on (here, a hypercube)
Definition bsplinebasis.hh:455
│ │ │ │ -
unsigned size() const
Number of shape functions in this finite element.
Definition bsplinebasis.hh:445
│ │ │ │ -
BSplineLocalCoefficients< dim > localCoefficients_
Definition bsplinebasis.hh:477
│ │ │ │ -
unsigned int size(int i) const
Number of degrees of freedom for one coordinate direction.
Definition bsplinebasis.hh:463
│ │ │ │ -
BSplineLocalBasis< GV, R > localBasis_
Definition bsplinebasis.hh:476
│ │ │ │ -
const BSplineLocalBasis< GV, R > & localBasis() const
Hand out a LocalBasis object.
Definition bsplinebasis.hh:427
│ │ │ │ -
Pre-basis for B-spline basis.
Definition bsplinebasis.hh:499
│ │ │ │ -
static constexpr size_type multiIndexBufferSize
Definition bsplinebasis.hh:567
│ │ │ │ -
std::array< unsigned, dim > elements_
Number of grid elements in the different coordinate directions.
Definition bsplinebasis.hh:1240
│ │ │ │ -
GridView gridView_
Definition bsplinebasis.hh:1242
│ │ │ │ -
double R
Definition bsplinebasis.hh:570
│ │ │ │ -
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:1085
│ │ │ │ -
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:781
│ │ │ │ -
std::array< unsigned int, dim > order_
Order of the B-spline for each space dimension.
Definition bsplinebasis.hh:1234
│ │ │ │ -
size_type dimension() const
Get the total dimension of the space spanned by this basis.
Definition bsplinebasis.hh:719
│ │ │ │ -
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:1128
│ │ │ │ -
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:1028
│ │ │ │ -
unsigned int size(size_t d) const
Number of shape functions in one direction.
Definition bsplinebasis.hh:774
│ │ │ │ -
GV GridView
The grid view that the FE space is defined on.
Definition bsplinebasis.hh:560
│ │ │ │ -
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:735
│ │ │ │ -
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:913
│ │ │ │ -
std::size_t size_type
Definition bsplinebasis.hh:561
│ │ │ │ -
static constexpr size_type maxMultiIndexSize
Definition bsplinebasis.hh:565
│ │ │ │ -
void update(const GridView &gv)
Update the stored grid view, to be called if the grid has changed.
Definition bsplinebasis.hh:686
│ │ │ │ -
void initializeIndices()
Initialize the global indices.
Definition bsplinebasis.hh:676
│ │ │ │ -
size_type size(const Dune::ReservedVector< ST, i > &prefix) const
Return number of possible values for next position in multi index.
Definition bsplinebasis.hh:712
│ │ │ │ -
unsigned int size() const
Total number of B-spline basis functions.
Definition bsplinebasis.hh:765
│ │ │ │ -
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition bsplinebasis.hh:680
│ │ │ │ -
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:1009
│ │ │ │ -
Node makeNode() const
Create tree node.
Definition bsplinebasis.hh:694
│ │ │ │ -
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:590
│ │ │ │ -
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:812
│ │ │ │ -
static constexpr size_type minMultiIndexSize
Definition bsplinebasis.hh:566
│ │ │ │ -
std::array< std::vector< double >, dim > knotVectors_
The knot vectors, one for each space dimension.
Definition bsplinebasis.hh:1237
│ │ │ │ -
size_type maxNodeSize() const
Get the maximal number of DOFs associated to node for any element.
Definition bsplinebasis.hh:725
│ │ │ │ -
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:642
│ │ │ │ -
LocalBasis class in the sense of dune-localfunctions, presenting the restriction of a B-spline patch ...
Definition bsplinebasis.hh:46
│ │ │ │ -
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:55
│ │ │ │ -
unsigned int order() const
Polynomial order of the shape functions.
Definition bsplinebasis.hh:140
│ │ │ │ -
std::size_t size() const
Return the number of basis functions on the current knot span.
Definition bsplinebasis.hh:147
│ │ │ │ -
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:97
│ │ │ │ -
void evaluateFunction(const FieldVector< D, dim > &in, std::vector< FieldVector< R, 1 > > &out) const
Evaluate all shape functions.
Definition bsplinebasis.hh:70
│ │ │ │ -
void evaluateJacobian(const FieldVector< D, dim > &in, std::vector< FieldMatrix< D, 1, dim > > &out) const
Evaluate Jacobian of all shape functions.
Definition bsplinebasis.hh:82
│ │ │ │ -
BSplineLocalBasis(const BSplinePreBasis< GV > &preBasis, const BSplineLocalFiniteElement< GV, R > &lFE)
Constructor with a given B-spline patch.
Definition bsplinebasis.hh:61
│ │ │ │ -
Attaches a shape function to an entity.
Definition bsplinebasis.hh:178
│ │ │ │ -
const LocalKey & localKey(std::size_t i) const
get i'th index
Definition bsplinebasis.hh:321
│ │ │ │ -
void init(const std::array< unsigned, dim > &sizes)
Definition bsplinebasis.hh:256
│ │ │ │ -
std::size_t size() const
number of coefficients
Definition bsplinebasis.hh:315
│ │ │ │ -
Local interpolation in the sense of dune-localfunctions, for the B-spline basis on tensor-product gri...
Definition bsplinebasis.hh:340
│ │ │ │ -
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition bsplinebasis.hh:344
│ │ │ │ -
Definition bsplinebasis.hh:1250
│ │ │ │ -
const FiniteElement & finiteElement() const
Return the LocalFiniteElement for the element we are bound to.
Definition bsplinebasis.hh:1274
│ │ │ │ -
typename GV::template Codim< 0 >::Entity Element
Definition bsplinebasis.hh:1256
│ │ │ │ -
const BSplinePreBasis< GV > * preBasis_
Definition bsplinebasis.hh:1290
│ │ │ │ -
Element element_
Definition bsplinebasis.hh:1293
│ │ │ │ -
void bind(const Element &e)
Bind to element.
Definition bsplinebasis.hh:1280
│ │ │ │ -
BSplineNode(const BSplinePreBasis< GV > *preBasis)
Definition bsplinebasis.hh:1259
│ │ │ │ -
const Element & element() const
Return current element, throw if unbound.
Definition bsplinebasis.hh:1265
│ │ │ │ -
FiniteElement finiteElement_
Definition bsplinebasis.hh:1292
│ │ │ │ -
std::size_t size_type
Definition bsplinebasis.hh:1255
│ │ │ │ -
Global basis for given pre-basis.
Definition defaultglobalbasis.hh:46
│ │ │ │ -
size_type size() const
Definition nodes.hh:142
│ │ │ │ -
void setSize(const size_type size)
Definition nodes.hh:164
│ │ │ │ -
│ │ │ │ +
An entity set for all entities of given codim in a grid view.
Definition gridviewentityset.hh:23
│ │ │ │ +
GridViewEntitySet(const GridView &gv)
Construct GridViewEntitySet for a GridView.
Definition gridviewentityset.hh:47
│ │ │ │ +
GV GridView
Definition gridviewentityset.hh:26
│ │ │ │ +
GridView::template Codim< codim >::Entity Element
Type of Elements contained in this EntitySet.
Definition gridviewentityset.hh:32
│ │ │ │ +
const_iterator end() const
Create an end iterator.
Definition gridviewentityset.hh:70
│ │ │ │ +
const GridView & gridView() const
Return the associated GridView.
Definition gridviewentityset.hh:76
│ │ │ │ +
Element value_type
Definition gridviewentityset.hh:38
│ │ │ │ +
const_iterator begin() const
Create a begin iterator.
Definition gridviewentityset.hh:64
│ │ │ │ +
GridView::template Codim< codim >::Iterator const_iterator
A forward iterator.
Definition gridviewentityset.hh:41
│ │ │ │ +
Element::Geometry::LocalCoordinate LocalCoordinate
Type of local coordinates with respect to the Element.
Definition gridviewentityset.hh:35
│ │ │ │ +
size_t size() const
Return number of Elements visited by an iterator.
Definition gridviewentityset.hh:58
│ │ │ │ +
Element::Geometry::GlobalCoordinate GlobalCoordinate
Definition gridviewentityset.hh:36
│ │ │ │ +
@ codim
Definition gridviewentityset.hh:28
│ │ │ │ +
bool contains(const Element &e) const
Return true if e is contained in the EntitySet.
Definition gridviewentityset.hh:52
│ │ │ │ +
const_iterator iterator
Same as const_iterator.
Definition gridviewentityset.hh:44
│ │ │ │