# -*- coding: utf-8 -*-
# noinspection PyUnresolvedReferences,PyRedeclaration
r"""
.. testsetup:: *

    None_or_custom_path_2 = './source/gallery/msepy_domains_and_meshes/msepy/backward_step2.png'
    None_or_custom_path_3 = './source/gallery/msepy_domains_and_meshes/msepy/backward_step3.png'

    import __init__ as ph
    from msepy.manifold.predefined.backward_step import _make_an_illustration
    _make_an_illustration(
        './source/gallery/msepy_domains_and_meshes/msepy/backward_step_illustration.png'
    )

.. testcleanup::

    pass


The backward step is a mesh (or domain) in :math:`\mathbb{R}^n`, :math:`n\in\left\lbrace2,3\right\rbrace`. The domain is
illustrated in the following figure.

.. figure:: backward_step_illustration.png
    :width: 100%

    The illustration of the backward step domain.

When :math:`n=3`, The domain is extended to the thrid axis, :math:`z`-axis, perpendicular to the plane. The parameters
are

.. autofunction:: msepy.manifold.predefined.backward_step.backward_step

Boundary units
==============

The domain is divided into three regions,

+---------------+--------------------------------------------------------------------------+
| region ``0``  |  bottom-right region, :math:`[x_1, (x_1+x_2)]\times[0, y_1]`             |
+---------------+--------------------------------------------------------------------------+
| region ``1``  |  top-right region , :math:`[x_1, (x_1+x_2)]\times[y_1, (y_1+y_2)]`       |
+---------------+--------------------------------------------------------------------------+
| region ``2``  |  top-left region, :math:`[0, x_1]\times[y_1, (y_1+y_2)]`                 |
+---------------+--------------------------------------------------------------------------+

Thus, for a 2-dimensional domain, it has 8 boundary units, i.e.,

>>> boundary_units_set = {
...     0: [1, 1, 1, 0],
...     1: [0, 1, 0, 1],
...     2: [1, 0, 1, 1],
... }

For example, for the bottom-right region ``0``, its left (:math:`x^-`), right (:math:`x^+`) and bottom
(:math:`y^-`) faces are boundary units, while its north (:math:`y^+`) face is not. So we have
``0: [1, 1, 1, 0]`` in the set.

And a 3-dimensional domain, it has 8 + 6
(:math:`2\times3` :math:`z`-direction) boundary units, i.e.,

>>> boundary_units_set = {
...     0: [1, 1, 1, 0, 1, 1],
...     1: [0, 1, 0, 1, 1, 1],
...     2: [1, 0, 1, 1, 1, 1],
... }

Examples
========

3d
--

Below codes will lead to a three-dimensional backward step mesh of :math:`3*5*5*5` elements (because the domain
is splited into 3 regions).

>>> ph.config.set_embedding_space_dim(3)
>>> manifold = ph.manifold(3)
>>> mesh = ph.mesh(manifold)
>>> msepy, obj = ph.fem.apply('msepy', locals())
>>> manifold = obj['manifold']
>>> mesh = obj['mesh']
>>> msepy.config(manifold)('backward_step')
>>> msepy.config(mesh)([5, 5, 5])
>>> mesh.visualize(saveto=None_or_custom_path_3)  # doctest: +ELLIPSIS
<Figure size ...

.. figure:: backward_step3.png
    :width: 60%

    The three-dimensional backward step mesh of :math:`3*5*5*5` elements.

2d
--

To make a two-dimensional backward step mesh of :math:`3*24*6` elements, just do

>>> ph.config.set_embedding_space_dim(2)
>>> manifold = ph.manifold(2)
>>> mesh = ph.mesh(manifold)
>>> msepy, obj = ph.fem.apply('msepy', locals())
>>> manifold = obj['manifold']
>>> mesh = obj['mesh']
>>> msepy.config(manifold)('backward_step')
>>> msepy.config(mesh)([24, 6])
>>> mesh.visualize(saveto=None_or_custom_path_2)  # doctest: +ELLIPSIS
<Figure size ...

.. figure:: backward_step2.png
    :width: 100%

    The two-dimensional backward step mesh of :math:`3*24*6` elements.


|

↩️  Back to :ref:`GALLERY-msepy-domains-and-meshes`.
"""


from msepy.manifold.predefined._helpers import _LinearTransformation

import matplotlib.pyplot as plt


def _make_an_illustration(saveto, x1=1, x2=1, y1=0.25, y2=0.25):
    """"""
    fig, ax = plt.subplots(figsize=(8, 6))
    ax.set_aspect('equal')
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    plt.tick_params(
        axis='both',
        which='both',
        left=False,
        bottom=False,
        labelbottom=False,
        labelleft=False
    )

    x = [x1, x1+x2, x1+x2, 0,     0, x1, x1]
    y = [0,  0,     y1+y2, y1+y2, y1, y1, 0]

    plt.plot(x, y, '-k', linewidth=0.8)
    plt.text(x1 + x2/2, 0.03, r'$x2$', ha='center')
    plt.text(x1/2, y1 + 0.03, r'$x1$', ha='center')
    plt.text(- 0.07, y1 + y2/2, r'$y2$', va='center')
    plt.text(x1 - 0.07, y1/2, r'$y1$', va='center')
    plt.text(x1, y1, r'$(x1, y1)$', va='bottom', ha='left')

    plt.xlabel(r"$x$")
    plt.ylabel(r"$y$")
    plt.savefig(saveto, bbox_inches='tight', dpi=200)
    plt.close()



[docs]
def backward_step(x1=1, x2=1, y1=0.25, y2=0.25, z=None, periodic=False):
    """

    Parameters
    ----------
    x1 : float, default=1
        See the illustration.
    x2 : float, default=1
        See the illustration.
    y1 : float, default=0.25
        See the illustration.
    y2 : float, default=0.25
        See the illustration.
    z : float, None, default=None
        When it is ``None``, it gives a two-dimensional domain. Otherwise, :math:`(z>0)`, it gives a three-dimensional
        one.
    periodic : bool, default=False
        When the domain is 3d, whether it is periodic along the ``z``-axis?

    """
    raise Exception(x1, x2, y1, y2, z, periodic)



def _backward_step(mf, x1=1, x2=1, y1=0.25, y2=0.25, z=None, x0=0, y0=0, periodic=False):
    r"""
    ^ y
    |
    |              x1                      x2
     __________________________________________________
    ||                    |                           |
    ||                    |                           |
    ||         r2         |           r1              |    y2
    ||                    |                           |
    ||____________________|___________________________|
    |            (x1, y1) |                           |
    |                     |                           |
    |                     |           r0              |    y1
    |                     |                           |
    | (x0,y0)             |___________________________|
    .--------------------------------------------------------------> x
    z

    Parameters
    ----------
    mf
    x1
    x2
    y1
    y2
    z

    Returns
    -------

    """
    assert mf.esd == mf.ndim, f"backward_step mesh only works for manifold.ndim == embedding space dimensions."
    assert mf.esd in (2, 3), f"backward_step mesh only works in 2-, 3-dimensions."
    esd = mf.esd
    if z is None:
        if esd == 2:
            z = 0
        else:
            z = 0.25
    else:
        pass
    if esd == 2:
        assert z == 0, f"for 2-d backward_step mesh, z must be 0."
    elif esd == 3:
        assert z > 0, f"for 3-d backward_step mesh, z must be greater than 0."
    else:
        raise NotImplementedError()

    if esd == 2:
        rm0 = _LinearTransformation(x0+x1, x0+x1+x2, y0+0,  y0+y1)
        rm1 = _LinearTransformation(x0+x1, x0+x1+x2, y0+y1, y0+y1+y2)
        rm2 = _LinearTransformation(x0+0,  x0+x1,    y0+y1, y0+y1+y2)
    elif esd == 3:
        rm0 = _LinearTransformation(x0+x1, x0+x1+x2, y0+0,  y0+y1,    0, z)
        rm1 = _LinearTransformation(x0+x1, x0+x1+x2, y0+y1, y0+y1+y2, 0, z)
        rm2 = _LinearTransformation(x0+0,  x0+x1,    y0+y1, y0+y1+y2, 0, z)
    else:
        raise Exception()

    if esd == 2:
        region_map = {
            0: [None, None, None, 1],
            1: [2,    None, 0,    None],
            2: [None, 1,    None, None],
        }
    elif esd == 3:
        if periodic:
            region_map = {
                0: [None, None, None, 1,    0, 0],
                1: [2,    None, 0,    None, 1, 1],
                2: [None, 1,    None, None, 2, 2],
            }
        else:
            region_map = {
                0: [None, None, None, 1,    None, None],
                1: [2,    None, 0,    None, None, None],
                2: [None, 1,    None, None, None, None],
            }
    else:
        raise Exception()

    mapping_dict = {
        0: rm0.mapping,  # region #0
        1: rm1.mapping,  # region #1
        2: rm2.mapping,  # region #2
    }

    Jacobian_matrix_dict = {
        0: rm0.Jacobian_matrix,  # region #1
        1: rm1.Jacobian_matrix,  # region #2
        2: rm2.Jacobian_matrix,  # region #3
    }

    if esd == 2:
        mtype_dict = {
            0: {'indicator': 'Linear', 'parameters': [f'x{x2}', f'y{y1}']},
            1: {'indicator': 'Linear', 'parameters': [f'x{x2}', f'y{y2}']},
            2: {'indicator': 'Linear', 'parameters': [f'x{x1}', f'y{y2}']},
        }
    elif esd == 3:
        mtype_dict = {
            0: {'indicator': 'Linear', 'parameters': [f'x{x2}', f'y{y1}', f'z{z}']},
            1: {'indicator': 'Linear', 'parameters': [f'x{x2}', f'y{y2}', f'z{z}']},
            2: {'indicator': 'Linear', 'parameters': [f'x{x1}', f'y{y2}', f'z{z}']}
        }
    else:
        raise Exception()

    return region_map, mapping_dict, Jacobian_matrix_dict, mtype_dict, None