writebasis.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
// SPDX-FileCopyrightText: Copyright © DUNE-FUFEM Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_WRITEBASIS_HH
#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_WRITEBASIS_HH
 
#include <type_traits>
#include <variant>
#include <string>
#include <algorithm>
 
#include <dune/common/exceptions.hh>
#include <dune/common/stringutility.hh>
#include <dune/common/version.hh>
 
#include <dune/grid/common/rangegenerators.hh>
 
#if DUNE_VERSION_GTE(DUNE_FUNCTIONS, 2, 11)
#include <dune/common/typetree/traversal.hh>
#else
#include <dune/typetree/traversal.hh>
#endif
 
 
#include <dune/functions/backends/istlvectorbackend.hh>
#include <dune/functions/backends/istlvectorfactory.hh>
#include <dune/functions/backends/containerfactory.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
 
#include <dune/vtk/vtkwriter.hh>
#include <dune/vtk/datacollectors/discontinuouslagrangedatacollector.hh>
 
#include <dune/fufem/backends/traversal.hh>
 
namespace Dune::Fufem {
 
  template<class Range=double, class Vector=std::monostate, class Basis>
  void writeBasis(const Basis& basis, std::string name, int order = -1)
  {
    // Create a vector suitable to hold coefficients for the basis
    auto vector = [&] {
      if constexpr (std::is_same_v<Vector, std::monostate>)
        return Dune::Functions::makeISTLVector<double>(basis.rootBasis().preBasis().containerDescriptor());
      else
      {
        auto v = Vector();
        Dune::Functions::istlVectorBackend(v).resize(basis);
        Dune::Functions::istlVectorBackend(v) = 0;
        return v;
      }
    }();
    auto vectorBackend = Dune::Functions::istlVectorBackend(vector);
 
    // Mark all DOFs of the subspace and deduce the appropriate order.
    // Since there's no mechanism for deducing a bitvector so far,
    // use a double vector to store flags.
    auto inSubspace = vector;
    auto inSubspaceBackend = Dune::Functions::istlVectorBackend(inSubspace);
    auto localView = basis.localView();
    for(const auto& element : Dune::elements(basis.gridView()))
    {
      localView.bind(element);
      Dune::TypeTree::forEachLeafNode(localView.tree(), [&](auto&& node, auto&& treePath) {
        for(auto i : Dune::range(node.size()))
          inSubspaceBackend[localView.index(node.localIndex(i))] = 1;
        if (order==-1)
          order = std::max(order, int(node.finiteElement().localBasis().order()));
      });
    }
 
    auto basisFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<double>(basis, vectorBackend);
    auto i_flat=0;
    Dune::Fufem::recursiveForEachVectorEntry(vector, [&](auto&& vector_i, const auto& i)
    {
      if (inSubspaceBackend[i]) {
        vectorBackend[i] = 1;
        auto vtkWriter = Dune::Vtk::UnstructuredGridWriter(Dune::Vtk::DiscontinuousLagrangeDataCollector(basis.gridView(), order));
        vtkWriter.addPointData(basisFunction, Dune::VTK::FieldInfo("basisfunction", Dune::VTK::FieldInfo::Type::scalar, 1));
        vtkWriter.write(Dune::formatString(name, i_flat));
        vectorBackend[i] = 0;
      }
      i_flat++;
    });
  }
 
  template<class Range=double, class Vector=std::monostate, class Basis, class Constraints>
  void writeBasis(const Basis& basis, std::string name, const Constraints& constraints, int order = -1)
  {
    // Create a vector suitable to hold coefficients for the basis
    auto vector = [&] {
      if constexpr (std::is_same_v<Vector, std::monostate>)
        return Dune::Functions::makeISTLVector<double>(basis.rootBasis().preBasis().containerDescriptor());
      else
      {
        auto v = Vector();
        Dune::Functions::istlVectorBackend(v).resize(basis);
        Dune::Functions::istlVectorBackend(v) = 0;
        return v;
      }
    }();
    auto vectorBackend = Dune::Functions::istlVectorBackend(vector);
 
    // Mark all DOFs of the subspace and deduce the appropriate order.
    // Since there's no mechanism for deducing a bitvector so far,
    // use a double vector to store flags.
    auto inSubspace = vector;
    auto inSubspaceBackend = Dune::Functions::istlVectorBackend(inSubspace);
    auto localView = basis.localView();
    for(const auto& element : Dune::elements(basis.gridView()))
    {
      localView.bind(element);
      Dune::TypeTree::forEachLeafNode(localView.tree(), [&](auto&& node, auto&& treePath) {
        for(auto i : Dune::range(node.size()))
          inSubspaceBackend[localView.index(node.localIndex(i))] = 1;
        if (order==-1)
          order = std::max(order, int(node.finiteElement().localBasis().order()));
      });
    }
 
    auto basisFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<double>(basis, vectorBackend);
    auto i_flat=0;
    Dune::Fufem::recursiveForEachVectorEntry(vector, [&](auto&& vector_i, const auto& i)
    {
      if (inSubspaceBackend[i]) {
        vectorBackend = 0;
        vectorBackend[i] = 1;
        constraints.interpolate(vector);
        auto vtkWriter = Dune::Vtk::UnstructuredGridWriter(Dune::Vtk::DiscontinuousLagrangeDataCollector(basis.gridView(), order));
        vtkWriter.addPointData(basisFunction, Dune::VTK::FieldInfo("basisfunction", Dune::VTK::FieldInfo::Type::scalar, 1));
        vtkWriter.write(Dune::formatString(name, i_flat));
      }
      i_flat++;
    });
  }
 
 
} // end namespace Dune::Fufem
 
#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_WRITEBASIS_HH