NSSC/Exercise_01/Solver.h

#pragma once

#include <functional>

#include "Matrix.h"

/** Stencil
 */
struct Stencil {
    const double C; /*< Center */
    const double B; /*< Bottom */
    const double L; /*< Left */
    const double U; /*< Up */
    const double R; /*< Right */
};

class Solver {
public:
    Solver(const size_t nx, const size_t ny,
        const double xmin, const double xmax,
        const double ymin, const double ymax,
        const double h, const double k,
        std::function<double(double, double)> fun,
        std::function<double(double)> gN,
        std::function<double(double)> gE,
        std::function<double(double)> gS,
        std::function<double(double)> gW
    ) :
        _iter{0},
        _nx{nx}, _ny{ny},
        _xmin{xmin}, _xmax{xmax},
        _ymin{ymin}, _ymax{ymax},
        _stencil{ 4.0 / (h * h) + 4.0 * k * k,  /* Center */
                 -1.0 / (h * h),                /* Bottom */
                 -1.0 / (h * h),                /* Left */
                 -1.0 / (h * h),                /* Up */
                 -1.0 / (h * h)},               /* Right */
        _rhs(nx, ny),
        _sol(nx, ny, 0.0),
        _tmp(nx, ny)
    {
        // Initialize Right Hand Size _rhs(x, y) = f(X(x), Y(y))
        // Note that the correspondence usual matrix indexing sceeme as
        // row/colums indices lead to an missplaces representation if the matrix
        // is printed in the usual format (x <-> rows, y <-> negative columns).
        // If this in entierly ignored and only considured in the case of
        // printing the matrix, everything is fine.
        for (size_t x = 0; x < nx; ++x) {
            for (size_t y = 0; y < ny; ++y) {
                _rhs(x, y) = fun(X(x), Y(y));
            }
        }

        // Set Dirichlet boundary conditions
        // North
        for (size_t x = 0; x < nx; ++x) { _sol(x, -1) = _tmp(x, -1) = gN(X(x)); }
        // East
        for (size_t y = 0; y < ny; ++y) { _sol(-1, y) = _tmp(-1, y) = gE(Y(y)); }
        // South
        for (size_t x = 0; x < nx; ++x) { _sol(x, 0) = _tmp(x, 0) = gS(X(x)); }
        // West
        for (size_t y = 0; y < ny; ++y) { _sol(0, y) = _tmp(0, y) = gW(Y(y)); }
    };

    /* Maps x/y index to x/y position in the [xmin, xmax] x [ymin, ymax] domain */
    double X(size_t x) const {
        double ratio = static_cast<double>(x) / static_cast<double>(_nx - 1);
        return _xmin + ratio * (_xmax - _xmin);
    }
    double Y(size_t y) const {
        double ratio = static_cast<double>(y) / static_cast<double>(_ny - 1);
        return _ymin + ratio * (_ymax - _ymin);
    }

    enum class Dir {
        North, East, South, West,
        N = North, E = East, S = South, W = West
    };

    MatrixView<double> boundary(enum Dir dir) {
        switch (dir) {
            case Dir::North:
                return Col<double>(_sol, -1);
            case Dir::East:
                return Row<double>(_sol, -1);
            case Dir::South:
                return Col<double>(_sol, 0);
            default:    // Dir::West
                return Row<double>(_sol, 0);
        }
    };

    /** Solution getter */
    Matrix<double>& solution() {
        return _sol;
    }

    /**
     * Performs a single Jacobian Iteration
     */
    void iterate() {
        double s = 1.0 / _stencil.C;

        for (size_t y = 1; y < _ny - 1; ++y) {
            for (size_t x = 1; x < _nx - 1; ++x) {
                _tmp(x, y) = s * (_rhs(x, y) - (_stencil.B * _sol(x, y - 1) +
                                                _stencil.L * _sol(x - 1, y) +
                                                _stencil.U * _sol(x, y + 1) +
                                                _stencil.R * _sol(x + 1, y)));
            }
        }

        std::swap(_tmp, _sol);
        _iter++;
    };

private:
    size_t _iter;           /*< Iteration count */
    const size_t _nx;       /*< Number of X-axis grid points */
    const size_t _ny;       /*< Number of Y-axis grid points */
    const double _xmin;     /*< Domain X-min (west border pos) */
    const double _xmax;     /*< Domain X-max (east border pos) */
    const double _ymin;     /*< Domain Y-min (south border pos) */
    const double _ymax;     /*< Domain Y-max (north border pos) */
    const Stencil _stencil; /*< Simple, + shaped stencil */
    Matrix<double> _rhs;    /*< Grid evaluated RHS of the PDE = f(x, y) */
    Matrix<double> _sol;    /*< Solution after _iter iterations */
    Matrix<double> _tmp;    /*< Temp. datablock, used in iterate() */
};