/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */
package edu.gcsc.vrl.modsim.pde;

import edu.gcsc.vrl.modsim.pde.predefined.NodeVisualization;

import eu.mihosoft.vrl.annotation.ComponentInfo;
import eu.mihosoft.vrl.annotation.MethodInfo;
import eu.mihosoft.vrl.annotation.ObjectInfo;
import eu.mihosoft.vrl.annotation.ParamInfo;
import eu.mihosoft.vrl.types.Function2D;

import java.io.Serializable;

/**
 * Simple "five point star" discretisation for the Poisson equation on the
 * unit square, discretised on an equidistant rectangular grid with N
 * rectangles in every spatial direction, so that the grid width is
 *
 * h := 1/N.
 *
 * Each grid point (x, y) has the representation x = ih, y = jh (0 <= i, j <= N).
 *
 * We traverse the grid in lexicographically "right-up" order
 * ie. one horizontal grid line from left to right after the other,
 * from bottom to top.
 *
 * NOTE: For the visualisation of the solution we have to assume this very same
 * ordering in 'LinearEquationSolver.java'!
 *
 * @author Ingo Heppner <Ingo.Heppner@gcsc.uni-frankfurt.de>
 */
@ComponentInfo(name = "AnisotropicFivePoint", category = "Tutorial/ModSim/PDE")
@ObjectInfo(name = "AnisotropicFivePoint")
public class AnisotropicFivePoint implements Serializable {

    private static final long serialVersionUID = 1;

    /**
     * computes the "five point stars"
     * @param f       2-D function describing the source term
     * @param phi     2-D function describing the Dirichlet boundary conditions
     * @param theGrid description of the underlying grid
     */
    @MethodInfo(valueName = "LES", valueOptions = "serialization=false")
    public LinearEquationSystem discretise(
            @ParamInfo(name = "f(x,y)") Function2D f,
            @ParamInfo(name = "phi(x,y)") Function2D phi,
            @ParamInfo(name = "Grid") Grid theGrid,
            @ParamInfo(name = "Label",
            options = "value=\"Five Point Star\"") String label) {

        if (theGrid.get_num_elems_per_dir() <= 0) {
            throw new IllegalArgumentException(
                    "FivePoint: N <= 0 is not allowed!");
        }

        // Notation after W. Hackbusch, "Iterative Loesung grosser schwachbesetzter
        // Gleichungssysteme", S. 18:
        // N: number of grid intervals in both spatial directions
        //    ==> number of grid points per direction, boundary inclusively, is N + 1,
        //        number of *inner* grid points per direction, without boundary, is N - 1.
        int Nx = theGrid.get_num_elems_in_x();
        int Ny = theGrid.get_num_elems_in_y();
        int NxPlus1 = Nx + 1;
        int NyPlus1 = Ny + 1;
        int nTotalPoints = theGrid.get_num_points();
        double hx = 1.0 / Nx;
        double hy = 1.0 / Ny;
        double hxSquared = hx * hx;
        double hySquared = hy * hy;

        int i;
        int j;
        int k = 0;

        /**
         * Latter:
         *
         *        [ (nTotalpoints-Nplus1)  ...      (nTotalPoint-1) ]                                ]
         *        [                        ...                      ]
         * y=j*h  [                        ...                      ]
         *        [ Nplus1 Nplus1+1        ...             2N  2N+1 ]
         *        [ k=0    1               ...        Nminus1     N ]
         *
         *             --> x=i*h
         *
         *
         */
        //int verbose = 8; // TMP
        // create equation system
        LinearEquationSystem les = new LinearEquationSystem(nTotalPoints);

        for (k = 0; k < NxPlus1; k++) {
            les.A[k][k] = 1;
            les.b[k] = phi.run(k * hx, 0.);
        }

        k = NxPlus1;

        for (j = 1; j < Ny; j++) {
            for (i = 0; i < NxPlus1; i++) {

                if (k % (NxPlus1) == 0 || k % (NxPlus1) == Nx) {
                    les.A[k][k] = 1;
                    les.b[k] = phi.run(i * hx, j * hy);
                } else {
                    les.A[k][k] = 2 * (hxSquared + hySquared) / (hxSquared * hySquared);
                    les.A[k][k - 1] = les.A[k][k + 1] = -1.0 / hxSquared;
                    les.A[k][k - NxPlus1] = les.A[k][k + NxPlus1] = -1.0 / hySquared;
                    les.b[k] = f.run(i * hx, j * hy);
                }
                k++;
            }
        }

        for (k = nTotalPoints - NxPlus1; k < nTotalPoints; k++) {
            i = k % NxPlus1;
            les.A[k][k] = 1;
            les.b[k] = phi.run(i * hx, Ny * hy);
        }


        // return computed linear equation system
        return les;
    }
}