assembler/interface.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers mueller@ibb.uni-stuttgart.de
// SPDX-License-Identifier: LGPL-3.0-or-later
 
#pragma once
#include <ranges>
 
#include <dune/common/referencehelper.hh>
 
#include <ikarus/assembler/dirichletbcenforcement.hh>
#include <ikarus/finiteelements/fehelper.hh>
#include <ikarus/finiteelements/ferequirements.hh>
#include <ikarus/utils/dirichletvalues.hh>
 
namespace Ikarus {
template <typename FEC, typename DV>
class FlatAssemblerBase
{
public:
  using FEContainerRaw = std::remove_cvref_t<FEC>; 
  using FERequirement  = typename FEContainerRaw::value_type::Requirement;
  using GlobalIndex = typename FEContainerRaw::value_type::GlobalIndex; 
  using Basis       = typename DV::Basis;                               
  using GridView    = typename Basis::GridView;                         
  using FEContainer = FEC;                                              
  using DirichletValuesType = DV;                          
  using SizeType = typename DirichletValuesType::SizeType; 
  using AffordanceCollectionType = AffordanceCollection<ScalarAffordance, VectorAffordance, MatrixAffordance>;
 
  template <typename FEContainer_ = FEContainer, typename DirichletValuesType_ = DirichletValuesType>
  FlatAssemblerBase(FEContainer_&& fes, DirichletValuesType_&& dirichletValues)
      : feContainer_{std::forward<FEContainer_>(fes)},
        dirichletValues_{std::forward<DirichletValuesType_>(dirichletValues)} {
    constraintsBelow_.reserve(dirichletValues_.size());
    size_t counter = 0;
    for (auto iv : std::ranges::iota_view{decltype(dirichletValues_.size())(0), dirichletValues_.size()}) {
      constraintsBelow_.emplace_back(counter);
      if (dirichletValues_.isConstrained(iv))
        ++counter;
    }
    fixedDofs_ = dirichletValues_.fixedDOFsize();
  }
 
  size_t reducedSize() { return dirichletValues_.size() - fixedDofs_; }
 
  size_t size() { return dirichletValues_.size(); }
 
  Eigen::VectorXd createFullVector(Eigen::Ref<const Eigen::VectorXd> reducedVector);
 
  const auto& finiteElements() const { return Dune::resolveRef(feContainer_); }
 
  auto& finiteElements() { return Dune::resolveRef(feContainer_); }
 
  const auto& dirichletValues() const { return Dune::resolveRef(dirichletValues_); }
 
  const auto& gridView() const { return Dune::resolveRef(dirichletValues_.basis().gridView()); }
 
  [[nodiscard]] size_t constraintsBelow(SizeType i) const { return constraintsBelow_[i]; }
 
  [[nodiscard]] bool isConstrained(SizeType i) const { return dirichletValues_.isConstrained(i); }
 
  [[nodiscard]] size_t estimateOfConnectivity() const {
    return dirichletValues_.basis().gridView().size(GridView::dimension) * 8;
  }
 
  void bind(const FERequirement& req, AffordanceCollectionType affordanceCollection,
            DBCOption dbcOption = DBCOption::Full) {
    req_         = std::make_optional<FERequirement>(req);
    affordances_ = std::make_optional<AffordanceCollectionType>(affordanceCollection);
    dBCOption_   = std::make_optional<DBCOption>(dbcOption);
  }
 
  void bind(const FERequirement& req) { req_ = std::make_optional<FERequirement>(req); }
 
  void bind(AffordanceCollectionType affordanceCollection) {
    affordances_ = std::make_optional<AffordanceCollectionType>(affordanceCollection);
  }
 
  void bind(DBCOption dbcOption) { dBCOption_ = std::make_optional<DBCOption>(dbcOption); }
 
  [[nodiscard]]
  bool bound() const {
    if (boundToRequirement() and boundToAffordanceCollection() and boundToDBCOption())
      return true;
    else
      DUNE_THROW(Dune::InvalidStateException, "The assembler is not bound to a requirement, affordance or dBCOption.");
  }
 
  [[nodiscard]]
  bool boundToRequirement() const {
    return req_.has_value();
  }
 
  [[nodiscard]]
  bool boundToAffordanceCollection() const {
    return affordances_.has_value();
  }
 
  [[nodiscard]]
  bool boundToDBCOption() const {
    return dBCOption_.has_value();
  }
 
  FERequirement& requirement() {
    if (req_.has_value())
      return req_.value();
    else
      DUNE_THROW(Dune::InvalidStateException, "The requirement can only be obtained after binding");
  }
 
  AffordanceCollectionType affordanceCollection() const {
    if (affordances_.has_value())
      return affordances_.value();
    else
      DUNE_THROW(Dune::InvalidStateException, "The affordance can only be obtained after binding");
  }
 
  DBCOption dBCOption() const {
    if (dBCOption_.has_value())
      return dBCOption_.value();
    else
      DUNE_THROW(Dune::InvalidStateException, "The dBCOption can only be obtained after binding");
  }
 
private:
  FEContainer feContainer_;
  DirichletValuesType dirichletValues_;
  std::optional<FERequirement> req_;
  std::optional<AffordanceCollectionType> affordances_;
  std::vector<size_t> constraintsBelow_{};
  size_t fixedDofs_{};
  std::optional<DBCOption> dBCOption_;
};
 
#ifndef DOXYGEN
template <class FEV, class DirichletValuesType>
FlatAssemblerBase(const FEV& fes,
                  const DirichletValuesType& dirichletValues) -> FlatAssemblerBase<FEV, DirichletValuesType>;
#endif
 
template <typename SA, typename FEC, typename DV, typename ST>
class ScalarAssembler
{
public:
  using ScalarAssemblerType = SA;
  using ScalarType          = ST;
  using DirichletValuesType = DV;
  using FEContainer         = FEC;
  using FEContainerRaw      = std::remove_cvref_t<FEC>; 
  using FERequirement       = typename FEContainerRaw::value_type::Requirement;
 
  ScalarType& scalar(const FERequirement& feRequirements, ScalarAffordance affordance) {
    return underlying().getScalarImpl(feRequirements, affordance);
  }
 
  ScalarType& scalar() {
    return underlying().getScalarImpl(underlying().requirement(),
                                      underlying().affordanceCollection().scalarAffordance());
  }
 
private:
  //> CRTP
  const auto& underlying() const { return static_cast<const ScalarAssemblerType&>(*this); }
  auto& underlying() { return static_cast<ScalarAssemblerType&>(*this); }
};
 
template <typename VA, typename FEC, typename DV, typename VT>
class VectorAssembler
{
public:
  using VectorAssemblerType = VA;
  using VectorType          = VT;
 
  using FEContainerRaw = std::remove_cvref_t<FEC>; 
  using FERequirement  = typename FEContainerRaw::value_type::Requirement;
  using GlobalIndex = typename FEContainerRaw::value_type::GlobalIndex; 
 
  using DirichletValuesType = DV;
  using FEContainer         = FEC;
 
  VectorType& vector(const FERequirement& feRequirements, VectorAffordance affordance,
                     DBCOption dbcOption = DBCOption::Full) {
    if (dbcOption == DBCOption::Raw) {
      return underlying().getRawVectorImpl(feRequirements, affordance);
    } else if (dbcOption == DBCOption::Reduced) {
      return underlying().getReducedVectorImpl(feRequirements, affordance);
    } else if (dbcOption == DBCOption::Full) {
      return underlying().getVectorImpl(feRequirements, affordance);
    }
    __builtin_unreachable();
  }
 
  VectorType& vector(DBCOption dbcOption) {
    return vector(underlying().requirement(), underlying().affordanceCollection().vectorAffordance(), dbcOption);
  }
 
  VectorType& vector() { return vector(underlying().dBCOption()); }
 
private:
  //> CRTP
  const auto& underlying() const { return static_cast<const VectorAssemblerType&>(*this); }
  auto& underlying() { return static_cast<VectorAssemblerType&>(*this); }
};
 
template <typename MA, typename FEC, typename DV, typename MT>
class MatrixAssembler
{
public:
  using MatrixAssemblerType = MA;
  using MatrixType          = MT;
 
  using FEContainerRaw = std::remove_cvref_t<FEC>; 
  using FERequirement  = typename FEContainerRaw::value_type::Requirement;
  using GlobalIndex = typename FEContainerRaw::value_type::GlobalIndex; 
 
  using DirichletValuesType = DV;
  using FEContainer         = FEC;
 
  MatrixType& matrix(const FERequirement& feRequirements, MatrixAffordance affordance,
                     DBCOption dbcOption = DBCOption::Full) {
    if (dbcOption == DBCOption::Raw) {
      return underlying().getRawMatrixImpl(feRequirements, affordance);
    } else if (dbcOption == DBCOption::Reduced) {
      return underlying().getReducedMatrixImpl(feRequirements, affordance);
    } else if (dbcOption == DBCOption::Full) {
      return underlying().getMatrixImpl(feRequirements, affordance);
    }
    __builtin_unreachable();
  }
 
  MatrixType& matrix(DBCOption dbcOption) {
    return matrix(underlying().requirement(), underlying().affordanceCollection().matrixAffordance(), dbcOption);
  }
 
  MatrixType& matrix() { return matrix(underlying().dBCOption()); }
 
private:
  //> CRTP
  const auto& underlying() const { return static_cast<const MatrixAssemblerType&>(*this); }
  auto& underlying() { return static_cast<MatrixAssemblerType&>(*this); }
};
} // namespace Ikarus