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 {
 
struct FemMatrix {
    const size_t rank;     
    const size_t size;     
    double* data;          
    const Solver* solver;  
 
    FemMatrix(const Solver* solver, size_t rank, size_t size)
        : rank(rank), size(size), data(aligned_malloc<double>(size)), solver(solver) {
        clear();
    }
 
    FemMatrix(const FemMatrix&) = delete;
 
    virtual ~FemMatrix() { aligned_free<double>(data); }
 
    virtual double& operator()(size_t r, size_t c) = 0;
 
    virtual void clear() {
        std::fill_n(data, size, 0.);
    }
 
    virtual void factorize() {}
 
    virtual void solverhs(DataVector<double>& B, DataVector<double>& X) = 0;
 
    void solve(DataVector<double>& B, DataVector<double>& X) {
        factorize();
        solverhs(B, X);
    }
 
    void solve(DataVector<double>& B) {
        solve(B, B);
    }
 
    virtual void mult(const DataVector<const double>& vector, DataVector<double>& result) = 0;
 
    virtual void addmult(const DataVector<const double>& vector, DataVector<double>& result) = 0;
 
    virtual void setBC(DataVector<double>& B, size_t r, double val) = 0;
 
    template <typename BoundaryConditonsT> void applyBC(const BoundaryConditonsT& bconds, DataVector<double>& 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);
    }
};
 
struct BandMatrix : FemMatrix {
    const size_t ld;       
    const size_t kd;       
 
    BandMatrix(const Solver* solver, size_t rank, size_t kd, size_t ld)
        : FemMatrix(solver, rank, rank * (ld + 1)), ld(ld), kd(kd) {}
 
    void setBC(DataVector<double>& B, size_t r, double val) override {
        B[r] = val;
        (*this)(r, r) = 1.;
        size_t start = (r > kd) ? r - kd : 0;
        size_t end = (r + kd < rank) ? r + kd + 1 : 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={}", rank, kd+1, size);
    }
};
 
 
}  // namespace plask
 
#endif