simpleassemblers.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2021-2024 The Ikarus Developers mueller@ibb.uni-stuttgart.de
// SPDX-License-Identifier: LGPL-3.0-or-later
 
#pragma once
 
#include <ranges>
#include <utility>
 
#include <dune/functions/backends/istlvectorbackend.hh>
 
#include <Eigen/Core>
#include <Eigen/Sparse>
 
#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 FEContainerType = std::conditional_t<std::is_reference_v<FEContainer>, const FEContainer, FEContainer>;
  using DirichletValuesType = DV; 
 
  FlatAssemblerBase(FEContainer&& fes, const DirichletValuesType& dirichletValues)
      : feContainer_{std::forward<FEContainer>(fes)},
        dirichletValues_{&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);
 
  auto& finiteElements() const { return feContainer_; }
 
  [[nodiscard]] size_t constraintsBelow(size_t i) const { return constraintsBelow_[i]; }
 
  [[nodiscard]] bool isConstrained(size_t i) const { return dirichletValues_->isConstrained(i); }
 
  [[nodiscard]] size_t estimateOfConnectivity() const {
    return dirichletValues_->basis().gridView().size(GridView::dimension) * 8;
  }
 
  using AffordanceCollectionType = AffordanceCollection<ScalarAffordance, VectorAffordance, MatrixAffordance>;
 
  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:
  FEContainerType feContainer_;
  const 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 T, class DirichletValuesType>
FlatAssemblerBase(T&& fes, const DirichletValuesType& dirichletValues) -> FlatAssemblerBase<T, DirichletValuesType>;
#endif
 
template <typename FEC, typename DV>
class ScalarAssembler : public FlatAssemblerBase<FEC, DV>
{
  using FEContainerRaw = std::remove_cvref_t<FEC>; 
  using Base           = FlatAssemblerBase<FEC, DV>;
 
public:
  using typename Base::Basis;
  using typename Base::DirichletValuesType;
  using typename Base::FEContainer;
  using typename Base::FERequirement;
  using typename Base::GlobalIndex;
 
  ScalarAssembler(FEContainer&& fes, const DirichletValuesType& dirichletValues)
      : FlatAssemblerBase<FEContainer, DirichletValuesType>(std::forward<FEContainer>(fes), dirichletValues) {}
 
  const double& scalar(const FERequirement& feRequirements, ScalarAffordance affordance) {
    return getScalarImpl(feRequirements, affordance);
  }
 
  const double& scalar() { return getScalarImpl(this->requirement(), this->affordanceCollection().scalarAffordance()); }
 
private:
  double& getScalarImpl(const FERequirement& feRequirements, ScalarAffordance affordance) {
    scal_ = 0.0;
    for (auto& fe : this->finiteElements()) {
      scal_ += calculateScalar(fe, feRequirements, affordance);
    }
    return scal_;
  }
 
  double scal_{0.0};
};
 
#ifndef DOXYGEN
template <class T, class DirichletValuesType>
ScalarAssembler(T&& fes, const DirichletValuesType& dirichletValues) -> ScalarAssembler<T, DirichletValuesType>;
#endif
 
template <typename FEC, typename DV>
class VectorFlatAssembler : public ScalarAssembler<FEC, DV>
{
  using FEContainerRaw = std::remove_cvref_t<FEC>; 
  using Base           = ScalarAssembler<FEC, DV>;
 
public:
  using typename Base::Basis;
  using typename Base::DirichletValuesType;
  using typename Base::FEContainer;
  using typename Base::FERequirement;
  using typename Base::GlobalIndex;
 
public:
  VectorFlatAssembler(FEContainer&& fes, const DirichletValuesType& dirichletValues)
      : ScalarAssembler<FEContainer, DirichletValuesType>(std::forward<FEContainer>(fes), dirichletValues) {}
 
  const Eigen::VectorXd& vector(const FERequirement& feRequirements, VectorAffordance affordance,
                                DBCOption dbcOption = DBCOption::Full) {
    if (dbcOption == DBCOption::Raw) {
      return getRawVectorImpl(feRequirements, affordance);
    } else if (dbcOption == DBCOption::Reduced) {
      return getReducedVectorImpl(feRequirements, affordance);
    } else if (dbcOption == DBCOption::Full) {
      return getVectorImpl(feRequirements, affordance);
    }
    __builtin_unreachable();
  }
 
  const Eigen::VectorXd& vector(DBCOption dbcOption) {
    return vector(this->requirement(), this->affordanceCollection().vectorAffordance(), dbcOption);
  }
 
  const Eigen::VectorXd& vector() { return vector(this->dBCOption()); }
 
private:
  void assembleRawVectorImpl(const FERequirement& feRequirements, VectorAffordance affordance,
                             Eigen::VectorXd& assemblyVec);
  Eigen::VectorXd& getRawVectorImpl(const FERequirement& feRequirements, VectorAffordance affordance);
  Eigen::VectorXd& getVectorImpl(const FERequirement& feRequirements, VectorAffordance affordance);
 
  Eigen::VectorXd& getReducedVectorImpl(const FERequirement& feRequirements, VectorAffordance affordance);
 
  Eigen::VectorXd vecRaw_{}; 
  Eigen::VectorXd vec_{};    
  Eigen::VectorXd vecRed_{}; 
};
 
#ifndef DOXYGEN
template <class T, class DirichletValuesType>
VectorFlatAssembler(T&& fes, const DirichletValuesType& dirichletValues) -> VectorFlatAssembler<T, DirichletValuesType>;
#endif
 
template <typename FEC, typename DV>
class SparseFlatAssembler : public VectorFlatAssembler<FEC, DV>
{
public:
  using FEContainerRaw = std::remove_cvref_t<FEC>; 
  using Base           = VectorFlatAssembler<FEC, DV>;
 
  using typename Base::Basis;
  using typename Base::DirichletValuesType;
  using typename Base::FEContainer;
  using typename Base::FERequirement;
  using typename Base::GlobalIndex;
 
  SparseFlatAssembler(FEContainer&& fes, const DirichletValuesType& dirichletValues)
      : VectorFlatAssembler<FEContainer, DirichletValuesType>(std::forward<FEContainer>(fes), dirichletValues) {}
 
  using GridView = typename Basis::GridView; 
 
  const Eigen::SparseMatrix<double>& matrix(const FERequirement& feRequirements, MatrixAffordance affordance,
                                            DBCOption dbcOption = DBCOption::Full) {
    if (dbcOption == DBCOption::Raw) {
      return getRawMatrixImpl(feRequirements, affordance);
    } else if (dbcOption == DBCOption::Reduced) {
      return getReducedMatrixImpl(feRequirements, affordance);
    } else if (dbcOption == DBCOption::Full) {
      return getMatrixImpl(feRequirements, affordance);
    }
    __builtin_unreachable();
  }
 
  const Eigen::SparseMatrix<double>& matrix(DBCOption dbcOption) {
    return matrix(this->requirement(), this->affordanceCollection().matrixAffordance(), dbcOption);
  }
 
  const Eigen::SparseMatrix<double>& matrix() { return matrix(this->dBCOption()); }
 
private:
  void assembleRawMatrixImpl(const FERequirement& feRequirements, MatrixAffordance affordance,
                             Eigen::SparseMatrix<double>& assemblyMat);
  Eigen::SparseMatrix<double>& getRawMatrixImpl(const FERequirement& feRequirements, MatrixAffordance affordance);
  Eigen::SparseMatrix<double>& getMatrixImpl(const FERequirement& feRequirements, MatrixAffordance affordance);
  Eigen::SparseMatrix<double>& getReducedMatrixImpl(const FERequirement& feRequirements, MatrixAffordance affordance);
 
  void createOccupationPattern(Eigen::SparseMatrix<double>& assemblyMat);
 
  void createReducedOccupationPattern(Eigen::SparseMatrix<double>& assemblyMat);
 
  void createLinearDOFsPerElement(Eigen::SparseMatrix<double>& assemblyMat);
 
  void createLinearDOFsPerElementReduced(Eigen::SparseMatrix<double>& assemblyMat);
 
  void preProcessSparseMatrix(Eigen::SparseMatrix<double>& assemblyMat);
 
  void preProcessSparseMatrixReduced(Eigen::SparseMatrix<double>& assemblyMat);
 
  Eigen::SparseMatrix<double> spMatRaw_;     
  Eigen::SparseMatrix<double> spMat_;        
  Eigen::SparseMatrix<double> spMatReduced_; 
  std::vector<std::vector<Eigen::Index>>
      elementLinearIndices_; 
  std::vector<std::vector<Eigen::Index>>
      elementLinearReducedIndices_; 
  bool sparsePreProcessorRaw_{}, sparsePreProcessor_{},
      sparsePreProcessorReduced_{}; 
};
 
#ifndef DOXYGEN
template <class FEC, class DV>
SparseFlatAssembler(FEC&& fes, const DV& dirichletValues) -> SparseFlatAssembler<FEC, DV>;
#endif
 
template <typename FEC, typename DV>
auto makeSparseFlatAssembler(FEC&& fes, const DV& dirichletValues) {
  return std::make_shared<SparseFlatAssembler<FEC, DV>>(std::forward<FEC>(fes), dirichletValues);
}
 
template <typename FEC, typename DV>
class DenseFlatAssembler : public VectorFlatAssembler<FEC, DV>
{
public:
  using FEContainerRaw = std::remove_cvref_t<FEC>;     
  using Base           = VectorFlatAssembler<FEC, DV>; 
 
  using typename Base::Basis;               
  using typename Base::DirichletValuesType; 
  using typename Base::FEContainer;         
  using typename Base::FERequirement;       
  using typename Base::GlobalIndex;         
 
  explicit DenseFlatAssembler(FEContainer&& fes, const DirichletValuesType& dirichletValues)
      : VectorFlatAssembler<FEContainer, DirichletValuesType>(std::forward<FEContainer>(fes), dirichletValues) {}
 
  const Eigen::MatrixXd& matrix(const FERequirement& feRequirements, MatrixAffordance affordance,
                                DBCOption dbcOption = DBCOption::Full) {
    if (dbcOption == DBCOption::Raw) {
      return getRawMatrixImpl(feRequirements, affordance);
    } else if (dbcOption == DBCOption::Reduced) {
      return getReducedMatrixImpl(feRequirements, affordance);
    } else if (dbcOption == DBCOption::Full) {
      return getMatrixImpl(feRequirements, affordance);
    }
    __builtin_unreachable();
  }
 
  const Eigen::MatrixXd& matrix(DBCOption dbcOption) {
    return matrix(this->requirement(), this->affordanceCollection().matrixAffordance(), dbcOption);
  }
 
  const Eigen::MatrixXd& matrix() { return matrix(this->dBCOption()); }
 
private:
  void assembleRawMatrixImpl(const FERequirement& feRequirements, MatrixAffordance affordance,
                             Eigen::MatrixXd& assemblyMat);
  Eigen::MatrixXd& getRawMatrixImpl(const FERequirement& feRequirements, MatrixAffordance affordance);
  Eigen::MatrixXd& getMatrixImpl(const FERequirement& feRequirements, MatrixAffordance affordance);
  Eigen::MatrixXd& getReducedMatrixImpl(const FERequirement& feRequirements, MatrixAffordance affordance);
 
  Eigen::MatrixXd matRaw_{}; 
  Eigen::MatrixXd mat_{};    
  Eigen::MatrixXd matRed_{}; 
};
 
#ifndef DOXYGEN
// https://en.cppreference.com/w/cpp/language/class_template_argument_deduction
template <class FEC, class DV>
DenseFlatAssembler(FEC&& fes, const DV& dirichletValues) -> DenseFlatAssembler<FEC, DV>;
#endif
 
template <typename FEC, typename DV>
auto makeDenseFlatAssembler(FEC&& fes, const DV& dirichletValues) {
  return std::make_shared<DenseFlatAssembler<FEC, DV>>(std::forward<FEC>(fes), dirichletValues);
}
} // namespace Ikarus
 
#include "simpleassemblers.inl"