ElementManager.h

Go to the documentation of this file.

#ifndef ELEMENT_MANAGER_H
#define ELEMENT_MANAGER_H
 
#include <algorithm>
#include <memory>
#include <string>
#include <vector>
 
#include "GenericElement.h"
#include "GenericMatrix.h"
#include "MatrixFactory.h"
#include "OneDimLinearElement.h"
 
template <typename T>
class ElementManager
{
    private:
        std::vector<std::unique_ptr<GenericElement<T>>> m_interior_cells{};
        std::vector<std::unique_ptr<GenericElement<T>>> m_boundary_cells{};
 
        void sortCells(std::vector<std::unique_ptr<GenericElement<T>>>& all_cells, int size);
        void flushMatrices();
 
    public:
        // It only makes sense for them to be public as I need to get references to these objects.
        std::shared_ptr<GenericMatrix<T>> m_capacitance_matrix{};
        std::shared_ptr<GenericMatrix<T>> m_stiffness_matrix{};
        std::valarray<T> m_forcing_vector{};
 
        ElementManager(int size, std::shared_ptr<MatrixFactory<T>> matrix_factory);
        ElementManager(int size, std::vector<std::unique_ptr<GenericElement<T>>>& all_cells,
                       std::shared_ptr<MatrixFactory<T>> matrix_factory);
 
        void assembleFemMatrices();
        void setCells(std::vector<std::unique_ptr<GenericElement<T>>>& all_cells);
};
 
template <typename T>
ElementManager<T>::ElementManager(const int size, std::shared_ptr<MatrixFactory<T>> matrix_factory)
    : m_capacitance_matrix(matrix_factory->createShared()),
      m_stiffness_matrix(matrix_factory->createShared())
{
    m_capacitance_matrix->buildAtRuntime(size, size);
    m_stiffness_matrix->buildAtRuntime(size, size);
    m_forcing_vector.resize(size);
}
 
template <typename T>
ElementManager<T>::ElementManager(const int size,
                                  std::vector<std::unique_ptr<GenericElement<T>>>& all_cells,
                                  std::shared_ptr<MatrixFactory<T>> matrix_factory)
{
    if (size == -1)
    {
        std::cerr << "[ElementManager]: Grid initialization failed, no matrix size could be retrived!\n";
        throw std::runtime_error(static_cast<std::string>("[ElementManager]: Grid initialization failed, no matrix size could be retrived!"));
    }
 
    m_capacitance_matrix = matrix_factory->createShared();
    m_stiffness_matrix = matrix_factory->createShared();
 
    // CSR build-up needs to be postponed until elements could be registered.
    if ((m_capacitance_matrix->getMatrixType() == MatrixType::Dense)
        && (m_stiffness_matrix->getMatrixType() == MatrixType::Dense))
    {
        m_capacitance_matrix->buildAtRuntime(size, size);
        m_stiffness_matrix->buildAtRuntime(size, size);
    }
    else if (m_capacitance_matrix->getMatrixType() == MatrixType::Dense
             || m_stiffness_matrix->getMatrixType() == MatrixType::Dense)
    {
        // This case should never be accessed, as both matrices are created together in
        std::cerr
            << "[ElementManager]: Non-supported mixing of dense and sparse matrix types!\nPleas use "
               "either just dense or just sparse matrices for the global matrices.\n";
        throw std::runtime_error(
            "[ElementManager]: Non-supported mixing of dense and sparse matrix types!\nPleas use "
            "either just dense or just sparse matrices for the global matrices.");
    }
    m_forcing_vector.resize(size);
    sortCells(all_cells, size);
}
 
template <typename T>
void ElementManager<T>::sortCells(std::vector<std::unique_ptr<GenericElement<T>>>& all_cells,
                                  int dim)
{
    // To build up a CSR matrix, the elements have to first be registered.
    if (m_capacitance_matrix->getMatrixType() == MatrixType::Csr
        && m_stiffness_matrix->getMatrixType() == MatrixType::Csr)
    {
        for (auto& cell_ptr : all_cells)
        {
            // All elements save the position to their entries in the matrix through reference
            // wrappers.
            cell_ptr->registerElementsCsr(*m_capacitance_matrix, *m_stiffness_matrix);
            if (cell_ptr->isBoundary())
            {
                for (auto& [key, boundary_ptr] : cell_ptr->getBoundaryCondition())
                {
                    boundary_ptr->registerElementsCsr(*m_capacitance_matrix, *m_stiffness_matrix);
                }
            }
        }
        m_capacitance_matrix->buildAtRuntime(dim, dim);
        m_stiffness_matrix->buildAtRuntime(dim, dim);
    }
    else if (m_capacitance_matrix->getMatrixType() == MatrixType::Csr
             || m_stiffness_matrix->getMatrixType() == MatrixType::Csr)
    {
        // This case should never be accessed, as both matrices are created together in
        std::cerr
            << "[ElementManager] Non-supported mixing of dense and sparse matrix types!\nPleas use "
               "either just dense or just sparse matrices for the global matrices.\n";
        throw std::runtime_error(
            "[ElementManager] Non-supported mixing of dense and sparse matrix types!\nPleas use "
            "either just dense or just sparse matrices for the global matrices.\n");
    }
    for (auto& cell_ptr : all_cells)
    {
        // All elements save the position to their entries in the matrix through reference wrappers.
        cell_ptr->getReferenceToElements(*m_capacitance_matrix, *m_stiffness_matrix,
                                         m_forcing_vector);
        if (cell_ptr->isBoundary())
        {
            for (auto& [key, boundary_ptr] : cell_ptr->getBoundaryCondition())
            {
                boundary_ptr->getReferenceToElements(*m_capacitance_matrix, *m_stiffness_matrix,
                                                     m_forcing_vector);
            }
            m_boundary_cells.push_back(std::move(cell_ptr));
            // The original cell list is made invalid here. This should also be reflected in Grid
            cell_ptr = nullptr;
        }
        else
        {
            m_interior_cells.push_back(std::move(cell_ptr));
            // The original cell list is made invalid here. This should also be reflected in Grid
            cell_ptr = nullptr;
        }
    }
}
 
template <typename T>
void ElementManager<T>::setCells(std::vector<std::unique_ptr<GenericElement<T>>>& all_cells)
{
    sortCells(all_cells);
}
 
template <typename T>
void ElementManager<T>::flushMatrices()
{
    m_stiffness_matrix->zeroMatrix();
    m_capacitance_matrix->zeroMatrix();
    m_forcing_vector *= 0;
}
 
template <typename T>
void ElementManager<T>::assembleFemMatrices()
{
    flushMatrices();
    for (auto& interior_cell_ptr : m_interior_cells)
    {
        interior_cell_ptr->writeMatrixVectorBlock();
    }
    for (auto& boundary_cell_ptr : m_boundary_cells)
    {
        boundary_cell_ptr->writeBoundaryCondition();
    }
}
 
#endif