matrix.hpp

Go to the documentation of this file.

/*
 * This file is part of PLaSK (https://plask.app) by Photonics Group at TUL
 * Copyright (c) 2023 Lodz University of Technology
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, version 3.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU General Public License for more details.
 */
#ifndef PLASK_COMMON_FEM_MATRIX_HPP
#define PLASK_COMMON_FEM_MATRIX_HPP
 
#include <plask/plask.hpp>
 
namespace plask {
 
template <typename T = double> struct FemMatrix {
    const size_t rank;     
    const size_t size;     
    T* data;               
    const Solver* solver;  
 
    FemMatrix(const Solver* solver, size_t rank, size_t size)
        : rank(rank), size(size), data(aligned_malloc<T>(size)), solver(solver) {
        clear();
    }
 
    FemMatrix(const FemMatrix&) = delete;
 
    virtual ~FemMatrix() { aligned_free<T>(data); }
 
    virtual T& operator()(size_t r, size_t c) = 0;
 
    virtual void clear() { std::fill_n(data, size, 0.); }
 
    virtual void factorize() {}
 
    virtual void solverhs(DataVector<T>& B, DataVector<T>& X) = 0;
 
    void solve(DataVector<T>& B, DataVector<T>& X) {
        factorize();
        solverhs(B, X);
    }
 
    void solve(DataVector<T>& B) { solve(B, B); }
 
    virtual void mult(const DataVector<const T>& vector, DataVector<T>& result) = 0;
 
    virtual void addmult(const DataVector<const T>& vector, DataVector<T>& result) = 0;
 
    virtual void setBC(DataVector<T>& B, size_t r, T val) = 0;
 
    template <typename BoundaryConditonsT> void applyBC(const BoundaryConditonsT& bconds, DataVector<T>& B) {
        // boundary conditions of the first kind
        for (auto cond : bconds) {
            for (auto r : cond.place) {
                setBC(B, r, cond.value);
            }
        }
    }
 
    virtual std::string describe() const { return format("rank={}, size={}", rank, size); }
};
 
template <typename T = double> struct BandMatrix : FemMatrix<T> {
    const size_t ld;  
    const size_t kd;  
 
    BandMatrix(const Solver* solver, size_t rank, size_t kd, size_t ld)
        : FemMatrix<T>(solver, rank, rank * (ld + 1)), ld(ld), kd(kd) {}
 
    void setBC(DataVector<T>& B, size_t r, T val) override {
        B[r] = val;
        (*this)(r, r) = 1.;
        size_t start = (r > kd) ? r - kd : 0;
        size_t end = (r + kd < this->rank) ? r + kd + 1 : this->rank;
        for (size_t c = start; c < r; ++c) {
            B[c] -= (*this)(r, c) * val;
            (*this)(r, c) = 0.;
        }
        for (size_t c = r + 1; c < end; ++c) {
            B[c] -= (*this)(r, c) * val;
            (*this)(r, c) = 0.;
        }
    }
 
    std::string describe() const override { return format("rank={}, bands={}, size={}", this->rank, kd + 1, this->size); }
};
 
}  // namespace plask
 
#endif