discreteglobalbasisfunction.hh 13.9 KB
Newer Older
1 2 3 4 5 6 7 8
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_FUNCTIONS_GRIDFUNCTIONS_DISCRETEGLOBALBASISFUNCTIONS_HH
#define DUNE_FUNCTIONS_GRIDFUNCTIONS_DISCRETEGLOBALBASISFUNCTIONS_HH

#include <memory>

#include <dune/common/shared_ptr.hh>
9
#include <dune/common/hybridutilities.hh>
10

11
#include <dune/functions/functionspacebases/hierarchicnodetorangemap.hh>
12
#include <dune/functions/functionspacebases/flatvectorview.hh>
13 14
#include <dune/functions/gridfunctions/gridviewentityset.hh>
#include <dune/functions/gridfunctions/gridfunction.hh>
15
#include <dune/functions/common/treedata.hh>
16 17
#include <dune/functions/backends/concepts.hh>
#include <dune/functions/backends/istlvectorbackend.hh>
18 19 20 21 22 23 24 25 26

namespace Dune {
namespace Functions {



/**
 * \brief A grid function induced by a global basis and a coefficient vector
 *
Carsten Gräser's avatar
Carsten Gräser committed
27 28
 * \ingroup FunctionImplementations
 *
29 30
 * This implements the grid function interface by combining a given global
 * basis and a coefficient vector. The part of the spanned space that should
31
 * be covered by the function is determined by a tree path that specifies the
32 33 34 35 36 37 38 39 40 41 42 43 44 45
 * corresponding local ansatz tree.
 *
 * This class supports mapping of subtrees to multi-component ranges,
 * vector-valued shape functions, and implicit product spaces given
 * by vector-valued coefficients. The mapping of these to the range
 * type is done via the following multistage procedure:
 *
 * 1.Each leaf node N in the local ansatz subtree is associated to an entry
 *   RE of the range-type via the given node-to-range-entry-map.
 *
 * Now let C be the coefficient block for a single basis function and
 * V the value of this basis function at the evaluation point. Notice
 * that both may be scalar, vector, matrix, or general container valued.
 *
46
 * 2.Each entry of C is associated with a flat index j via flatVectorView.
47 48
 *   This is normally a lexicographic index. The total scalar dimension according
 *   to those flat indices is dim(C).
49
 * 3.Each entry of V is associated with a flat index k via flatVectorView.
50 51
 *   This is normally a lexicographic index. The total scalar dimension according
 *   to those flat indices dim(V).
52
 * 4.Each entry of RE is associated with a flat index k via flatVectorView.
53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69
 *   This is normally a lexicographic index. The total scalar dimension according
 *   to those flat indices dim(RE).
 * 5.Via those flat indices we now interpret C,V, and RE as vectors and compute the diadic
 *   product (C x V). The entries of this product are mapped to the flat indices for
 *   RE lexicographically. I.e. we set
 *
 *     RE[j*dim(V)+k] = C[j] * V[k]
 *
 * Hence the range entry RE must have dim(RE) = dim(C)*dim(V).
 *
 * \tparam B Type of lobal basis
 * \tparam TP Type of tree path specifying the requested subtree of ansatz functions
 * \tparam V Type of coefficient vectors
 * \tparam NTRE Type of node-to-range-entry-map that associates each leaf node in the local ansatz subtree with an entry in the range type
 * \tparam R Range type of this function
 */
template<typename B, typename TP, typename V,
70
  typename NTRE = HierarchicNodeToRangeMap,
71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94
  typename R = typename V::value_type>
class DiscreteGlobalBasisFunction
{
public:
  using Basis = B;
  using TreePath = TP;
  using Vector = V;

  using GridView = typename Basis::GridView;
  using EntitySet = GridViewEntitySet<GridView, 0>;
  using Tree = typename Basis::LocalView::Tree;
  using SubTree = typename TypeTree::ChildForTreePath<Tree, TreePath>;
  using NodeToRangeEntry = NTRE;

  using Domain = typename EntitySet::GlobalCoordinate;
  using Range = R;

  using LocalDomain = typename EntitySet::LocalCoordinate;
  using Element = typename EntitySet::Element;

  using Traits = Imp::GridFunctionTraits<Range(Domain), EntitySet, DefaultDerivativeTraits, 16>;

  class LocalFunction
  {
95
    using LocalView = typename Basis::LocalView;
96 97
    using size_type = typename SubTree::size_type;

98 99 100 101 102 103 104 105
    template<class Node>
    using LocalBasisRange = typename Node::FiniteElement::Traits::LocalBasisType::Traits::RangeType;

    template<class Node>
    using NodeData = typename std::vector<LocalBasisRange<Node>>;

    using ShapeFunctionValueContainer = TreeData<SubTree, NodeData, true>;

106 107
    struct LocalEvaluateVisitor
      : public TypeTree::TreeVisitor
108
      , public TypeTree::StaticTraversal
109
    {
110
      LocalEvaluateVisitor(const LocalDomain& x, Range& y, const LocalView& localView, const Vector& coefficients, const NodeToRangeEntry& nodeToRangeEntry, ShapeFunctionValueContainer& shapeFunctionValueContainer):
111 112
        x_(x),
        y_(y),
113
        localView_(localView),
114
        coefficients_(coefficients),
115 116
        nodeToRangeEntry_(nodeToRangeEntry),
        shapeFunctionValueContainer_(shapeFunctionValueContainer)
117 118 119 120 121
      {}

      template<typename Node, typename TreePath>
      void leaf(Node& node, TreePath treePath)
      {
122 123
        auto&& fe = node.finiteElement();
        auto&& localBasis = fe.localBasis();
124

125
        auto&& shapeFunctionValues = shapeFunctionValueContainer_[node];
126 127 128
        localBasis.evaluateFunction(x_, shapeFunctionValues);

        // Get range entry associated to this node
129
        auto re = flatVectorView(nodeToRangeEntry_(node, treePath, y_));
130 131 132 133


        for (size_type i = 0; i < localBasis.size(); ++i)
        {
134
          auto&& multiIndex = localView_.index(node.localIndex(i));
135 136

          // Get coefficient associated to i-th shape function
137
          auto c = flatVectorView(coefficients_[multiIndex]);
138 139

          // Get value of i-th shape function
140 141 142
          auto v = flatVectorView(shapeFunctionValues[i]);


143 144

          // Notice that the range entry re, the coefficient c, and the shape functions
145
          // value v may all be scalar, vector, matrix, or general container valued.
146 147
          // The matching of their entries is done via the multistage procedure described
          // in the class documentation of DiscreteGlobalBasisFunction.
148 149 150
          auto&& dimC = c.size();
          auto dimV = v.size();
          assert(dimC*dimV == re.size());
151 152
          for(size_type j=0; j<dimC; ++j)
          {
153
            auto&& c_j = c[j];
154
            for(size_type k=0; k<dimV; ++k)
155
              re[j*dimV + k] += c_j*v[k];
156 157 158 159 160 161
          }
        }
      }

      const LocalDomain& x_;
      Range& y_;
162
      const LocalView& localView_;
163 164
      const Vector& coefficients_;
      const NodeToRangeEntry& nodeToRangeEntry_;
165
      ShapeFunctionValueContainer& shapeFunctionValueContainer_;
166 167 168 169 170 171 172 173 174 175 176 177
    };


  public:

    using GlobalFunction = DiscreteGlobalBasisFunction;
    using Domain = LocalDomain;
    using Range = GlobalFunction::Range;
    using Element = GlobalFunction::Element;

    LocalFunction(const DiscreteGlobalBasisFunction& globalFunction)
      : globalFunction_(&globalFunction)
178
      , localView_(globalFunction.basis().localView())
179
    {
180 181
      // Here we assume that the tree can be accessed, traversed,
      // and queried for tree indices even in unbound state.
182
      subTree_ = &TypeTree::child(localView_.tree(), globalFunction_->treePath());
183
      shapeFunctionValueContainer_.init(*subTree_);
184
//      localDoFs_.reserve(localView_.maxSize());
185 186
    }

187 188
    LocalFunction(const LocalFunction& other)
      : globalFunction_(other.globalFunction_)
189
      , localView_(other.localView_)
190
      , bound_(other.bound_)
191 192 193
    {
      // Here we assume that the tree can be accessed, traversed,
      // and queried for tree indices even in unbound state.
194
      subTree_ = &TypeTree::child(localView_.tree(), globalFunction_->treePath());
195 196 197 198 199 200
      shapeFunctionValueContainer_.init(*subTree_);
    }

    LocalFunction operator=(const LocalFunction& other)
    {
      globalFunction_ = other.globalFunction_;
201 202
      localView_ = other.localView_;
      subTree_ = &TypeTree::child(localView_.tree(), globalFunction_->treePath());
203 204 205 206 207 208

      // Here we assume that the tree can be accessed, traversed,
      // and queried for tree indices even in unbound state.
      shapeFunctionValueContainer_.init(*subTree_);
    }

209 210 211 212 213 214 215 216
    /**
     * \brief Bind LocalFunction to grid element.
     *
     * You must call this method before evaluate()
     * and after changes to the coefficient vector.
     */
    void bind(const Element& element)
    {
217
      localView_.bind(element);
218
      bound_ = true;
219 220 221 222

      // Read dofs associated to bound element
//      localDoFs_.resize(subTree.size());
//      for (size_type i = 0; i < subTree.size(); ++i)
223
//        localDoFs_[i] = globalFunction_->dofs()[localView.index(i)];
224 225 226 227
    }

    void unbind()
    {
228
      localView_.unbind();
229 230 231 232 233 234 235 236
      bound_ = false;
    }

    /**
     * \brief Check if LocalFunction is already bound to an element.
     */
    bool bound() const {
      return bound_;
237 238 239 240 241 242 243 244 245 246 247 248 249 250 251
    }

    /**
     * \brief Evaluate LocalFunction at bound element.
     *
     * The result of this method is undefined if you did
     * not call bind() beforehand or changed the coefficient
     * vector after the last call to bind(). In the latter case
     * you have to call bind() again in order to make operator()
     * usable.
     */
    Range operator()(const Domain& x) const
    {
      auto y = Range(0);

252
      LocalEvaluateVisitor localEvaluateVisitor(x, y, localView_, globalFunction_->dofs(), globalFunction_->nodeToRangeEntry(), shapeFunctionValueContainer_);
253 254 255 256 257 258 259
      TypeTree::applyToTree(*subTree_, localEvaluateVisitor);

      return y;
    }

    const Element& localContext() const
    {
260
      return localView_.element();
261 262 263 264 265 266 267 268 269 270
    }

    friend typename Traits::LocalFunctionTraits::DerivativeInterface derivative(const LocalFunction& t)
    {
      DUNE_THROW(NotImplemented,"not implemented");
    }

  private:

    const DiscreteGlobalBasisFunction* globalFunction_;
271
    LocalView localView_;
272 273

    mutable ShapeFunctionValueContainer shapeFunctionValueContainer_;
274 275
//    std::vector<typename V::value_type> localDoFs_;
    const SubTree* subTree_;
276 277

    bool bound_ = false;
278 279
  };

280 281
  template<class B_T, class V_T, class NTRE_T>
  DiscreteGlobalBasisFunction(B_T && basis, const TreePath& treePath, V_T && coefficients, NTRE_T&& nodeToRangeEntry) :
282
    entitySet_(basis.gridView()),
283
    basis_(wrap_or_move(std::forward<B_T>(basis))),
284
    treePath_(treePath),
285 286
    coefficients_(wrap_or_move(std::forward<V_T>(coefficients))),
    nodeToRangeEntry_(wrap_or_move(std::forward<NTRE_T>(nodeToRangeEntry)))
287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327
  {}

  DiscreteGlobalBasisFunction(std::shared_ptr<const Basis> basis, const TreePath& treePath, std::shared_ptr<const V> coefficients, std::shared_ptr<const NodeToRangeEntry> nodeToRangeEntry) :
    entitySet_(basis->gridView()),
    basis_(basis),
    treePath_(treePath),
    coefficients_(coefficients),
    nodeToRangeEntry_(nodeToRangeEntry)
  {}

  const Basis& basis() const
  {
    return *basis_;
  }

  const TreePath& treePath() const
  {
    return treePath_;
  }

  const V& dofs() const
  {
    return *coefficients_;
  }

  const NodeToRangeEntry& nodeToRangeEntry() const
  {
    return *nodeToRangeEntry_;
  }

  // TODO: Implement this using hierarchic search
  Range operator() (const Domain& x) const
  {
    DUNE_THROW(NotImplemented,"not implemented");
  }

  friend typename Traits::DerivativeInterface derivative(const DiscreteGlobalBasisFunction& t)
  {
    DUNE_THROW(NotImplemented,"not implemented");
  }

328 329 330 331 332 333 334 335 336 337 338 339 340 341 342
  /**
   * \brief Construct local function from a DiscreteGlobalBasisFunction
   *
   * \ingroup FunctionImplementations
   *
   * The obtained local function satisfies the concept
   * `Dune::Functions::Concept::LocalFunction` and must be bound
   * to an entity from the entity set of the DiscreteGlobalBasisFunction
   * before it can be used.
   *
   * Notice that the local function stores a reference to the
   * global DiscreteGlobalBasisFunction. Hence calling any method
   * of the local function after the DiscreteGlobalBasisFunction
   * exceeded its life time leads to undefined behavior.
   */
343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363
  friend LocalFunction localFunction(const DiscreteGlobalBasisFunction& t)
  {
    return LocalFunction(t);
  }

  /**
   * \brief Get associated EntitySet
   */
  const EntitySet& entitySet() const
  {
    return entitySet_;
  }

private:

  EntitySet entitySet_;
  std::shared_ptr<const Basis> basis_;
  const TreePath treePath_;
  std::shared_ptr<const V> coefficients_;
  std::shared_ptr<const NodeToRangeEntry> nodeToRangeEntry_;
};
364 365 366



367 368 369 370 371 372 373 374 375 376
/**
 * \brief Construction of local functions from a temporary DiscreteGlobalBasisFunction (forbidden)
 *
 * Since a DiscreteGlobalBasisFunction::LocalFunction stores a reference
 * to the global DiscreteGlobalBasisFunction its life time is bound to
 * the latter. Hence construction from a temporary DiscreteGlobalBasisFunction
 * would lead to a dangling reference and is thus forbidden/deleted.
 *
 * \ingroup FunctionImplementations
 */
377 378
template<typename... TT>
void localFunction(DiscreteGlobalBasisFunction<TT...>&& t) = delete;
379 380 381 382



template<typename R, typename B, typename TP, typename V>
383
auto makeDiscreteGlobalBasisFunction(B&& basis, const TP& treePath, V&& vector)
384
{
385
  using Basis = std::decay_t<B>;
386
  using NTREM = HierarchicNodeToRangeMap;
387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403

  // Small helper functions to wrap vectors using istlVectorBackend
  // if they do not already satisfy the VectorBackend interface.
  auto toConstVectorBackend = [&](auto&& v) -> decltype(auto) {
    return Dune::Hybrid::ifElse(models<Concept::ConstVectorBackend<Basis>, decltype(v)>(),
    [&](auto id) -> decltype(auto) {
      return std::forward<decltype(v)>(v);
    }, [&](auto id) -> decltype(auto) {
      return istlVectorBackend(v);
    });
  };

  using Vector = std::decay_t<decltype(toConstVectorBackend(std::forward<V>(vector)))>;
  return DiscreteGlobalBasisFunction<Basis, TP, Vector, NTREM, R>(
      std::forward<B>(basis),
      treePath,
      toConstVectorBackend(std::forward<V>(vector)),
404
      HierarchicNodeToRangeMap());
405 406 407 408
}



409 410 411 412 413 414 415 416
template<typename R, typename B, typename V>
auto makeDiscreteGlobalBasisFunction(B&& basis, V&& vector)
{
  return makeDiscreteGlobalBasisFunction<R>(std::forward<B>(basis), TypeTree::hybridTreePath(), std::forward<V>(vector));
}



417 418 419 420
} // namespace Functions
} // namespace Dune

#endif // DUNE_FUNCTIONS_GRIDFUNCTIONS_DISCRETEGLOBALBASISFUNCTIONS_HH