/*
 * Copyright (c) NM LTD.
 * https://nm.dev/
 * 
 * THIS SOFTWARE IS LICENSED, NOT SOLD.
 * 
 * YOU MAY USE THIS SOFTWARE ONLY AS DESCRIBED IN THE LICENSE.
 * IF YOU ARE NOT AWARE OF AND/OR DO NOT AGREE TO THE TERMS OF THE LICENSE,
 * DO NOT USE THIS SOFTWARE.
 * 
 * THE SOFTWARE IS PROVIDED "AS IS", WITH NO WARRANTY WHATSOEVER,
 * EITHER EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION,
 * ANY WARRANTIES OF ACCURACY, ACCESSIBILITY, COMPLETENESS,
 * FITNESS FOR A PARTICULAR PURPOSE, MERCHANTABILITY, NON-INFRINGEMENT, 
 * TITLE AND USEFULNESS.
 * 
 * IN NO EVENT AND UNDER NO LEGAL THEORY,
 * WHETHER IN ACTION, CONTRACT, NEGLIGENCE, TORT, OR OTHERWISE,
 * SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR
 * ANY CLAIMS, DAMAGES OR OTHER LIABILITIES,
 * ARISING AS A RESULT OF USING OR OTHER DEALINGS IN THE SOFTWARE.
 */
package dev.nm.nmj;

import dev.nm.algebra.linear.vector.doubles.Vector;
import dev.nm.analysis.differentialequation.pde.finitedifference.PDESolutionGrid2D;
import dev.nm.analysis.differentialequation.pde.finitedifference.PDESolutionTimeSpaceGrid1D;
import dev.nm.analysis.differentialequation.pde.finitedifference.PDESolutionTimeSpaceGrid2D;
import dev.nm.analysis.differentialequation.pde.finitedifference.elliptic.dim2.IterativeCentralDifference;
import dev.nm.analysis.differentialequation.pde.finitedifference.elliptic.dim2.PoissonEquation2D;
import dev.nm.analysis.differentialequation.pde.finitedifference.hyperbolic.dim1.ExplicitCentralDifference1D;
import dev.nm.analysis.differentialequation.pde.finitedifference.hyperbolic.dim1.WaveEquation1D;
import dev.nm.analysis.differentialequation.pde.finitedifference.hyperbolic.dim2.ExplicitCentralDifference2D;
import dev.nm.analysis.differentialequation.pde.finitedifference.hyperbolic.dim2.WaveEquation2D;
import dev.nm.analysis.differentialequation.pde.finitedifference.parabolic.dim1.heatequation.CrankNicolsonHeatEquation1D;
import dev.nm.analysis.differentialequation.pde.finitedifference.parabolic.dim1.heatequation.HeatEquation1D;
import dev.nm.analysis.differentialequation.pde.finitedifference.parabolic.dim2.AlternatingDirectionImplicitMethod;
import dev.nm.analysis.differentialequation.pde.finitedifference.parabolic.dim2.HeatEquation2D;
import dev.nm.analysis.function.rn2r1.AbstractBivariateRealFunction;
import dev.nm.analysis.function.rn2r1.AbstractTrivariateRealFunction;
import dev.nm.analysis.function.rn2r1.BivariateRealFunction;
import dev.nm.analysis.function.rn2r1.TrivariateRealFunction;
import dev.nm.analysis.function.rn2r1.univariate.AbstractUnivariateRealFunction;
import dev.nm.number.DoubleUtils;
import static java.lang.Math.*;

/**
 * Numerical Methods Using Java: For Data Science, Analysis, and Engineering
 *
 * @author haksunli
 * @see
 * https://www.amazon.com/Numerical-Methods-Using-Java-Engineering/dp/1484267966
 * https://nm.dev/
 */
public class Chapter8 {

    public static void main(String[] args) {
        System.out.println("Chapter 8 demos");

        Chapter8 chapter8 = new Chapter8();
        chapter8.solve_Poisson_equation_0010();
        chapter8.solve_Poisson_equation_0020();
        chapter8.solve_wave_equation_1D();
        chapter8.solve_wave_equation_2D();
        chapter8.solve_heat_equation_1D();
        chapter8.solve_heat_equation_2D();
    }

    private void solve_heat_equation_2D() {
        System.out.println("solve a 2-dimensional heat equation");

        // solution domain
        final double a = 4.0, b = 4.0;
        // time domain
        final double T = 5000;

        // heat equation coefficient
        final double beta = 1e-4;

        // initial condition
        final BivariateRealFunction f
                = new AbstractBivariateRealFunction() {
            @Override
            public double evaluate(double x1, double x2) {
                return 0.;
            }
        };

        // boundary condition
        final TrivariateRealFunction g
                = new AbstractTrivariateRealFunction() {
            @Override
            public double evaluate(double t, double x, double y) {
                return exp(y) * cos(x) - exp(x) * cos(y);
            }
        };
        final HeatEquation2D PDE = new HeatEquation2D(beta, T, a, b, f, g);

        AlternatingDirectionImplicitMethod adi
                = new AlternatingDirectionImplicitMethod(1e-5);
        PDESolutionTimeSpaceGrid2D soln = adi.solve(PDE, 50, 39, 39);

        int t = 50;
        int x = 1;
        int y = 1;
        System.out.println(String.format("u(%d,%d,%d) = %f", t, x, y, soln.u(t, x, y)));
        t = 50;
        x = 1;
        y = 16;
        System.out.println(String.format("u(%d,%d,%d) = %f", t, x, y, soln.u(t, x, y)));
        t = 50;
        x = 1;
        y = 31;
        System.out.println(String.format("u(%d,%d,%d) = %f", t, x, y, soln.u(t, x, y)));

        t = 50;
        x = 16;
        y = 1;
        System.out.println(String.format("u(%d,%d,%d) = %f", t, x, y, soln.u(t, x, y)));
        t = 50;
        x = 16;
        y = 16;
        System.out.println(String.format("u(%d,%d,%d) = %f", t, x, y, soln.u(t, x, y)));
        t = 50;
        x = 16;
        y = 31;
        System.out.println(String.format("u(%d,%d,%d) = %f", t, x, y, soln.u(t, x, y)));

        t = 50;
        x = 31;
        y = 1;
        System.out.println(String.format("u(%d,%d,%d) = %f", t, x, y, soln.u(t, x, y)));
        t = 50;
        x = 31;
        y = 16;
        System.out.println(String.format("u(%d,%d,%d) = %f", t, x, y, soln.u(t, x, y)));
        t = 50;
        x = 31;
        y = 31;
        System.out.println(String.format("u(%d,%d,%d) = %f", t, x, y, soln.u(t, x, y)));
    }

    private void solve_heat_equation_1D() {
        System.out.println("solve a 1-dimensional heat equation");

        HeatEquation1D pde = new HeatEquation1D(
                1e-5, // heat equation coefficient
                1., 6000., // solution domain bounds
                new AbstractUnivariateRealFunction() {
            @Override
            public double evaluate(double x) {
                return 2.0 * x + sin(2.0 * PI * x); // initial condition
            }
        },
                0., new AbstractUnivariateRealFunction() {
            @Override
            public double evaluate(double t) {
                return 0.; // boundary condition at x = 0
            }
        },
                0., new AbstractUnivariateRealFunction() {
            @Override
            public double evaluate(double t) {
                return 2.; // boundary condition at x = 6000
            }
        });

        // c_k are 0 for Dirichlet boundary conditions
        int m = 50;
        int n = 39;

        PDESolutionTimeSpaceGrid1D soln
                = new CrankNicolsonHeatEquation1D().solve(pde, m, n);

        int t = 0;
        int x = 1;
        System.out.println(String.format("u(%d,%d) = %f", t, x, soln.u(t, x)));
        t = 0;
        x = 16;
        System.out.println(String.format("u(%d,%d) = %f", t, x, soln.u(t, x)));
        t = 0;
        x = 31;
        System.out.println(String.format("u(%d,%d) = %f", t, x, soln.u(t, x)));

        t = 15;
        x = 1;
        System.out.println(String.format("u(%d,%d) = %f", t, x, soln.u(t, x)));
        t = 15;
        x = 16;
        System.out.println(String.format("u(%d,%d) = %f", t, x, soln.u(t, x)));
        t = 15;
        x = 31;
        System.out.println(String.format("u(%d,%d) = %f", t, x, soln.u(t, x)));

        t = 30;
        x = 1;
        System.out.println(String.format("u(%d,%d) = %f", t, x, soln.u(t, x)));
        t = 30;
        x = 16;
        System.out.println(String.format("u(%d,%d) = %f", t, x, soln.u(t, x)));
        t = 30;
        x = 31;
        System.out.println(String.format("u(%d,%d) = %f", t, x, soln.u(t, x)));

        t = 45;
        x = 1;
        System.out.println(String.format("u(%d,%d) = %f", t, x, soln.u(t, x)));
        t = 45;
        x = 16;
        System.out.println(String.format("u(%d,%d) = %f", t, x, soln.u(t, x)));
        t = 45;
        x = 31;
        System.out.println(String.format("u(%d,%d) = %f", t, x, soln.u(t, x)));
    }

    private void solve_wave_equation_2D() {
        System.out.println("solve a 2-dimensional wave equation");

        double c2 = 1. / 4; // wave speed squared
        double T = 2., a = 2., b = 2.; // the solution domain bounds
        WaveEquation2D pde = new WaveEquation2D(
                c2, T, a, b,
                new AbstractBivariateRealFunction() {
            @Override
            public double evaluate(double x, double y) {
                return 0.1 * sin(PI * x) * sin(PI * y / 2.);
            }
        },
                new AbstractBivariateRealFunction() {
            @Override
            public double evaluate(double x, double y) {
                return 0.;
            }
        });

        int m = 40; // dt = T/m
        int n = 39; // dx = a/n
        int p = 39; // dy = b/p
        PDESolutionTimeSpaceGrid2D soln = new ExplicitCentralDifference2D().solve(pde, m, n, p);

        int t = 40; // t index
        int x = 1; // x index
        int y = 1; // y index
        System.out.println(String.format("u(%d,%d,%d) = %f", t, x, y, soln.u(t, x, y)));
        t = 40; // t index
        x = 1; // x index
        y = 16; // y index
        System.out.println(String.format("u(%d,%d,%d) = %f", t, x, y, soln.u(t, x, y)));
        t = 40; // t index
        x = 1; // x index
        y = 31; // y index
        System.out.println(String.format("u(%d,%d,%d) = %f", t, x, y, soln.u(t, x, y)));

        t = 40; // t index
        x = 16; // x index
        y = 1; // y index
        System.out.println(String.format("u(%d,%d,%d) = %f", t, x, y, soln.u(t, x, y)));
        t = 40; // t index
        x = 16; // x index
        y = 16; // y index
        System.out.println(String.format("u(%d,%d,%d) = %f", t, x, y, soln.u(t, x, y)));
        t = 40; // t index
        x = 16; // x index
        y = 31; // y index
        System.out.println(String.format("u(%d,%d,%d) = %f", t, x, y, soln.u(t, x, y)));

        t = 40; // t index
        x = 31; // x index
        y = 1; // y index
        System.out.println(String.format("u(%d,%d,%d) = %f", t, x, y, soln.u(t, x, y)));
        t = 40; // t index
        x = 31; // x index
        y = 16; // y index
        System.out.println(String.format("u(%d,%d,%d) = %f", t, x, y, soln.u(t, x, y)));
        t = 40; // t index
        x = 31; // x index
        y = 31; // y index
        System.out.println(String.format("u(%d,%d,%d) = %f", t, x, y, soln.u(t, x, y)));
    }

    private void solve_wave_equation_1D() {
        System.out.println("solve a 1-dimensional wave equation");

        final double c2 = 4.0; // c^2
        final double T = 1.0; // time upper bond
        final double a = 2.0; // x upper bound
        WaveEquation1D pde
                = new WaveEquation1D(
                        c2,
                        T,
                        a,
                        new AbstractUnivariateRealFunction() {
                    @Override
                    public double evaluate(double x) {
                        return 0.1 * sin(PI * x); // 0.1 * sin(π x)
                    }
                },
                        new AbstractUnivariateRealFunction() {
                    @Override
                    public double evaluate(double x) {
                        return 0.2 * PI * sin(PI * x); // 0.2π * sin(π x)
                    }
                });

        int m = 80; // dt = T/m
        int n = 39; // dx = a/n
        PDESolutionTimeSpaceGrid1D soln = new ExplicitCentralDifference1D().solve(pde, m, n);

        int t = 0; // time index
        int x = 1; // x index
        System.out.println(String.format("u(%d,%d) = %f", t, x, soln.u(t, x)));
        t = 0;
        x = 16;
        System.out.println(String.format("u(%d,%d) = %f", t, x, soln.u(t, x)));
        t = 0;
        x = 31;
        System.out.println(String.format("u(%d,%d) = %f", t, x, soln.u(t, x)));

        t = 20;
        x = 1;
        System.out.println(String.format("u(%d,%d) = %f", t, x, soln.u(t, x)));
        t = 20;
        x = 16;
        System.out.println(String.format("u(%d,%d) = %f", t, x, soln.u(t, x)));
        t = 20;
        x = 31;
        System.out.println(String.format("u(%d,%d) = %f", t, x, soln.u(t, x)));

        t = 40;
        x = 1;
        System.out.println(String.format("u(%d,%d) = %f", t, x, soln.u(t, x)));
        t = 40;
        x = 16;
        System.out.println(String.format("u(%d,%d) = %f", t, x, soln.u(t, x)));
        t = 40;
        x = 31;
        System.out.println(String.format("u(%d,%d) = %f", t, x, soln.u(t, x)));

        t = 60;
        x = 1;
        System.out.println(String.format("u(%d,%d) = %f", t, x, soln.u(t, x)));
        t = 60;
        x = 16;
        System.out.println(String.format("u(%d,%d) = %f", t, x, soln.u(t, x)));
        t = 60;
        x = 31;
        System.out.println(String.format("u(%d,%d) = %f", t, x, soln.u(t, x)));
    }

    private void solve_Poisson_equation_0020() {
        System.out.println("solve a 2-dimensional Poisson equation");

        BivariateRealFunction ZERO // a constant zero function, f = 0
                = new AbstractBivariateRealFunction() {
            @Override
            public double evaluate(double x, double y) {
                return 0;
            }

            @Override
            public Double evaluate(Vector x) {
                return 0.;
            }
        };

        // the boundary conditions
        final double EPSION = 1e-8;
        BivariateRealFunction g = new AbstractBivariateRealFunction() {
            @Override
            public double evaluate(double x, double y) {
                if (DoubleUtils.isZero(x, EPSION) || DoubleUtils.isZero(y, EPSION)) {
                    return 0;
                } else if (DoubleUtils.equal(x, 0.5, EPSION)) {
                    return 200. * y;
                } else if (DoubleUtils.equal(y, 0.5, EPSION)) {
                    return 200. * x;
                }

                // not reachable; don't matter
                return Double.NaN;
            }
        };

        double a = 0.5; // width of the x-dimension
        double b = 0.5; // height of the y-dimension
        PoissonEquation2D pde = new PoissonEquation2D(a, b, ZERO, g);
        IterativeCentralDifference solver = new IterativeCentralDifference(
                EPSION, // precision
                40); // max number of iterations
        PDESolutionGrid2D soln = solver.solve(pde, 4, 4);
        int k = 1, j = 1; // node indices
        double u_11 = soln.u(k, j);
        System.out.println(String.format("u_%d,%d = u(%f,%f): %f", k, j, soln.x(k), soln.y(j), u_11));
        k = 1;
        j = 2;
        double u_12 = soln.u(k, j);
        System.out.println(String.format("u_%d,%d = u(%f,%f): %f", k, j, soln.x(k), soln.y(j), u_12));
        k = 2;
        j = 1;
        double u_21 = soln.u(k, j);
        System.out.println(String.format("u_%d,%d = u(%f,%f): %f", k, j, soln.x(k), soln.y(j), u_21));
        k = 2;
        j = 2;
        double u_22 = soln.u(k, j);
        System.out.println(String.format("u_%d,%d = u(%f,%f): %f", k, j, soln.x(k), soln.y(j), u_22));
        k = 3;
        j = 3;
        double u_33 = soln.u(k, j);
        System.out.println(String.format("u_%d,%d = u(%f,%f): %f", k, j, soln.x(k), soln.y(j), u_33));
        k = 4;
        j = 4;
        double u_44 = soln.u(k, j);
        System.out.println(String.format("u_%d,%d = u(%f,%f): %f", k, j, soln.x(k), soln.y(j), u_44));
        k = 5;
        j = 5;
        double u_55 = soln.u(k, j);
        System.out.println(String.format("u_%d,%d = u(%f,%f): %f", k, j, soln.x(k), soln.y(j), u_55));
    }

    private void solve_Poisson_equation_0010() {
        System.out.println("solve a 2-dimensional Poisson equation");

        BivariateRealFunction ZERO // a constant zero function, f = 0
                = new AbstractBivariateRealFunction() {
            @Override
            public double evaluate(double x, double y) {
                return 0;
            }

            @Override
            public Double evaluate(Vector x) {
                return 0.;
            }
        };

        // the boundary conditions
        BivariateRealFunction g = new AbstractBivariateRealFunction() {
            @Override
            public double evaluate(double x, double y) {
                return log((1. + x) * (1. + x) + y * y);
            }
        };

        double a = 1.; // width of the x-dimension
        double b = 1.; // height of the y-dimension
        PoissonEquation2D pde = new PoissonEquation2D(a, b, ZERO, g);
        IterativeCentralDifference solver = new IterativeCentralDifference(
                1e-8, // precision
                40); // max number of iterations
        PDESolutionGrid2D soln = solver.solve(pde, 2, 2);
        int k = 1, j = 1; // node indices
        double u_11 = soln.u(k, j); // x = 0.3, y = 0.3
        System.out.println(String.format("u_%d,%d = u(%f,%f): %f", k, j, soln.x(k), soln.y(j), u_11));
        k = 1;
        j = 2;
        double u_12 = soln.u(k, j); // x = 0.3, y = 0.6
        System.out.println(String.format("u_%d,%d = u(%f,%f): %f", k, j, soln.x(k), soln.y(j), u_12));
        k = 2;
        j = 1;
        double u_21 = soln.u(k, j); // x = 0.6, y = 0.3
        System.out.println(String.format("u_%d,%d = u(%f,%f): %f", k, j, soln.x(k), soln.y(j), u_21));
        k = 2;
        j = 2;
        double u_22 = soln.u(k, j); // x = 0.6, y = 0.6
        System.out.println(String.format("u_%d,%d = u(%f,%f): %f", k, j, soln.x(k), soln.y(j), u_22));
    }

}