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/utils/dirichletvalues.hh>
 
namespace Ikarus {
  template <typename FEContainer_, typename DirichletValuesType_>
  class FlatAssemblerBase {
  public:
    using FEContainerRaw    = std::remove_cvref_t<FEContainer_>;  
    using FERequirementType = typename FEContainerRaw::value_type::FERequirementType;
    using GlobalIndex     = typename FEContainerRaw::value_type::GlobalIndex;  
    using Basis           = typename DirichletValuesType_::Basis;              
    using GridView        = typename Basis::GridView;                          
    using FEContainer     = FEContainer_;  
    using FEContainerType = std::conditional_t<std::is_reference_v<FEContainer>, const FEContainer, FEContainer>;
    using DirichletValuesType = DirichletValuesType_;  
 
    FlatAssemblerBase(FEContainer&& fes, const DirichletValuesType& p_dirichletValues)
        : feContainer{std::forward<FEContainer>(fes)}, dirichletValues{&p_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;
    DirichletValuesType const* 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 FEContainer_, typename DirichletValuesType_>
  class ScalarAssembler : public FlatAssemblerBase<FEContainer_, DirichletValuesType_> {
    using FEContainerRaw = std::remove_cvref_t<FEContainer_>;  
    using Base           = FlatAssemblerBase<FEContainer_, DirichletValuesType_>;
 
  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 += fe.calculateScalar(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 FEContainer_, typename DirichletValuesType_>
  class VectorFlatAssembler : public ScalarAssembler<FEContainer_, DirichletValuesType_> {
    using FEContainerRaw = std::remove_cvref_t<FEContainer_>;  
    using Base           = ScalarAssembler<FEContainer_, DirichletValuesType_>;
 
  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 FEContainer_, typename DirichletValuesType_>
  class SparseFlatAssembler : public VectorFlatAssembler<FEContainer_, DirichletValuesType_> {
  public:
    using FEContainerRaw = std::remove_cvref_t<FEContainer_>;  
    using Base           = VectorFlatAssembler<FEContainer_, DirichletValuesType_>;
 
    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 T, class DirichletValuesType>
  SparseFlatAssembler(T&& fes, const DirichletValuesType& dirichletValues_)
      -> SparseFlatAssembler<T, DirichletValuesType>;
#endif
 
  template <typename FEContainer_, typename DirichletValuesType_>
  class DenseFlatAssembler : public VectorFlatAssembler<FEContainer_, DirichletValuesType_> {
  public:
    using FEContainerRaw = std::remove_cvref_t<FEContainer_>;  
    using Base           = VectorFlatAssembler<FEContainer_, DirichletValuesType_>;  
 
    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 T, class DirichletValuesType>
  DenseFlatAssembler(T&& fes, const DirichletValuesType& dirichletValues_)
      -> DenseFlatAssembler<T, DirichletValuesType>;
#endif
}  // namespace Ikarus
 
#include "simpleassemblers.inl"