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 <mutex>
#include <ranges>
#include <utility>
 
#include <dune/functions/backends/istlvectorbackend.hh>
 
#include <Eigen/Core>
#include <Eigen/Sparse>
 
#include <ikarus/finiteelements/fehelper.hh>
#include <ikarus/utils/dirichletvalues.hh>
 
namespace Ikarus {
template <typename FEC, typename DV>
class FlatAssemblerBase
{
public:
  using FEContainerRaw    = std::remove_cvref_t<FEC>; 
  using FERequirementType = typename FEContainerRaw::value_type::FERequirementType;
  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;
  }
 
private:
  FEContainerType feContainer_;
  const DirichletValuesType* dirichletValues_;
  std::vector<size_t> constraintsBelow_{};
  size_t fixedDofs_{};
};
 
#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::FERequirementType;
  using typename Base::GlobalIndex;
 
  ScalarAssembler(FEContainer&& fes, const DirichletValuesType& dirichletValues)
      : FlatAssemblerBase<FEContainer, DirichletValuesType>(std::forward<FEContainer>(fes), dirichletValues) {}
 
  const double& getScalar(const FERequirementType& feRequirements) { return getScalarImpl(feRequirements); }
 
private:
  double& getScalarImpl(const FERequirementType& feRequirements) {
    scal_ = 0.0;
    for (auto& fe : this->finiteElements()) {
      scal_ += calculateScalar(fe, feRequirements);
    }
    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::FERequirementType;
  using typename Base::GlobalIndex;
 
public:
  VectorFlatAssembler(FEContainer&& fes, const DirichletValuesType& dirichletValues)
      : ScalarAssembler<FEContainer, DirichletValuesType>(std::forward<FEContainer>(fes), dirichletValues) {}
 
  const Eigen::VectorXd& getRawVector(const FERequirementType& feRequirements) {
    return getRawVectorImpl(feRequirements);
  }
 
  const Eigen::VectorXd& getVector(const FERequirementType& feRequirements) { return getVectorImpl(feRequirements); }
 
  const Eigen::VectorXd& getReducedVector(const FERequirementType& feRequirements) {
    return getReducedVectorImpl(feRequirements);
  }
 
private:
  void assembleRawVectorImpl(const FERequirementType& feRequirements, Eigen::VectorXd& assemblyVec);
  Eigen::VectorXd& getRawVectorImpl(const FERequirementType& feRequirements);
  Eigen::VectorXd& getVectorImpl(const FERequirementType& feRequirements);
 
  Eigen::VectorXd& getReducedVectorImpl(const FERequirementType& feRequirements);
 
  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::FERequirementType;
  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>& getRawMatrix(const FERequirementType& feRequirements) {
    return getRawMatrixImpl(feRequirements);
  }
 
  const Eigen::SparseMatrix<double>& getMatrix(const FERequirementType& feRequirements) {
    return getMatrixImpl(feRequirements);
  }
 
  const Eigen::SparseMatrix<double>& getReducedMatrix(const FERequirementType& feRequirements) {
    return getReducedMatrixImpl(feRequirements);
  }
 
private:
  void assembleRawMatrixImpl(const FERequirementType& feRequirements, Eigen::SparseMatrix<double>& assemblyMat);
  Eigen::SparseMatrix<double>& getRawMatrixImpl(const FERequirementType& feRequirements);
  Eigen::SparseMatrix<double>& getMatrixImpl(const FERequirementType& feRequirements);
  Eigen::SparseMatrix<double>& getReducedMatrixImpl(const FERequirementType& feRequirements);
 
  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_; 
  std::once_flag 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>
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::FERequirementType;   
  using typename Base::GlobalIndex;         
 
  explicit DenseFlatAssembler(FEContainer&& fes, const DirichletValuesType& dirichletValues)
      : VectorFlatAssembler<FEContainer, DirichletValuesType>(std::forward<FEContainer>(fes), dirichletValues) {}
 
  const Eigen::MatrixXd& getRawMatrix(const FERequirementType& feRequirements) {
    return getRawMatrixImpl(feRequirements);
  }
 
  const Eigen::MatrixXd& getMatrix(const FERequirementType& feRequirements) { return getMatrixImpl(feRequirements); }
 
  const Eigen::MatrixXd& getReducedMatrix(const FERequirementType& feRequirements) {
    return getReducedMatrixImpl(feRequirements);
  }
 
private:
  void assembleRawMatrixImpl(const FERequirementType& feRequirements, Eigen::MatrixXd& assemblyMat);
  Eigen::MatrixXd& getRawMatrixImpl(const FERequirementType& feRequirements);
  Eigen::MatrixXd& getMatrixImpl(const FERequirementType& feRequirements);
  Eigen::MatrixXd& getReducedMatrixImpl(const FERequirementType& feRequirements);
 
  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
} // namespace Ikarus
 
#include "simpleassemblers.inl"