import math
import os
import re
import time
from pathlib import Path
from subprocess import run

import gmsh
import matplotlib.pyplot as plt
import numpy as np
import ogstools as ot
import pyvista as pv

E = 1.0
nu = 0.15
Gc = 1.0
h = 0.01
a0 = 0.05  # half of the initial crack length

phasefield_model = "AT1"  # AT1/AT2

# file's name
prj_name = "Kregime_Propagating.prj"
meshname = "mesh_full_pf"

out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
if not out_dir.exists():
    out_dir.mkdir(parents=True)

def mesh_generation(lc, lc_fine):
    """
    lc ... characteristic length for coarse meshing
    lc_fine ... characteristic length for fine meshing
    """
    L = 4.0  # Length
    H = 4.0  # Height
    b = 0.4  # Length/Height of subdomain with fine mesh

    # Before using any functions in the Python API, Gmsh must be initialized
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 1)
    gmsh.model.add("rectangle")

    # Dimensions
    dim0 = 0
    dim1 = 1
    dim2 = 2

    # Points
    gmsh.model.geo.addPoint(-L / 2, -H / 2, 0, lc, 1)
    gmsh.model.geo.addPoint(L / 2, -H / 2, 0, lc, 2)
    gmsh.model.geo.addPoint(L / 2, H / 2, 0, lc, 3)
    gmsh.model.geo.addPoint(-L / 2, H / 2, 0, lc, 4)
    gmsh.model.geo.addPoint(-b, -b - lc_fine / 2, 0, lc_fine, 5)
    gmsh.model.geo.addPoint(b, -b - lc_fine / 2, 0, lc_fine, 6)
    gmsh.model.geo.addPoint(b, b + lc_fine / 2, 0, lc_fine, 7)
    gmsh.model.geo.addPoint(-b, b + lc_fine / 2, 0, lc_fine, 8)

    # Lines
    gmsh.model.geo.addLine(1, 2, 1)
    gmsh.model.geo.addLine(2, 3, 2)
    gmsh.model.geo.addLine(3, 4, 3)
    gmsh.model.geo.addLine(4, 1, 4)
    gmsh.model.geo.addLine(5, 6, 5)
    gmsh.model.geo.addLine(6, 7, 6)
    gmsh.model.geo.addLine(7, 8, 7)
    gmsh.model.geo.addLine(8, 5, 8)

    # Line loops
    gmsh.model.geo.addCurveLoop([1, 2, 3, 4], 1)
    gmsh.model.geo.addCurveLoop([5, 6, 7, 8], 2)

    # Add plane surfaces defined by one or more curve loops.
    gmsh.model.geo.addPlaneSurface([1, 2], 1)
    gmsh.model.geo.addPlaneSurface([2], 2)

    gmsh.model.geo.synchronize()

    # Prepare structured grid
    gmsh.model.geo.mesh.setTransfiniteCurve(
        6, math.ceil(2 * b / lc_fine + 2), "Progression", 1
    )
    gmsh.model.geo.mesh.setTransfiniteCurve(
        8, math.ceil(2 * b / lc_fine + 2), "Progression", 1
    )
    gmsh.model.geo.mesh.setTransfiniteSurface(2, "Alternate")

    gmsh.model.geo.mesh.setRecombine(dim2, 1)
    gmsh.model.geo.mesh.setRecombine(dim2, 2)

    gmsh.model.geo.synchronize()

    # Physical groups
    P1 = gmsh.model.addPhysicalGroup(dim0, [1])
    gmsh.model.setPhysicalName(dim0, P1, "P1")

    P2 = gmsh.model.addPhysicalGroup(dim0, [2])
    gmsh.model.setPhysicalName(dim0, P2, "P2")

    Bottom = gmsh.model.addPhysicalGroup(dim1, [1])
    gmsh.model.setPhysicalName(dim1, Bottom, "Bottom")

    Right = gmsh.model.addPhysicalGroup(dim1, [2])
    gmsh.model.setPhysicalName(dim1, Right, "Right")

    Top = gmsh.model.addPhysicalGroup(dim1, [3])
    gmsh.model.setPhysicalName(dim1, Top, "Top")

    Left = gmsh.model.addPhysicalGroup(dim1, [4])
    gmsh.model.setPhysicalName(dim1, Left, "Left")

    Computational_domain = gmsh.model.addPhysicalGroup(dim2, [1, 2])
    gmsh.model.setPhysicalName(dim2, Computational_domain, "Computational_domain")
    gmsh.model.geo.synchronize()

    output_file = f"{out_dir}/" + meshname + ".msh"
    gmsh.model.mesh.generate(dim2)
    gmsh.write(output_file)
    gmsh.finalize()

def pre_processing(h, a0):
    mesh = pv.read(f"{out_dir}/mesh_full_pf_domain.vtu")
    phase_field = np.ones((len(mesh.points), 1))
    pv.set_plot_theme("document")

    for node_id, x in enumerate(mesh.points):
        if (
            (mesh.center[0] - x[0]) <= a0 + 0.001 * h
            and (mesh.center[0] - x[0]) >= -a0 - 0.001 * h
            and (mesh.center[1] - x[1]) < h / 2 + 0.001 * h
            and (mesh.center[1] - x[1]) > -h / 2 - 0.001 * h
        ):
            phase_field[node_id] = 0.0

    mesh.point_data["phase-field"] = phase_field
    mesh.save(f"{out_dir}/mesh_full_pf_OGS_pf_ic.vtu")

pv.set_plot_theme("document")
pv.set_jupyter_backend("static")


def Hydraulic_Fracturing_Toughness_Dominated_numerical(h, phasefield_model):
    # mesh properties
    ls = 2 * h
    # generate prefix from properties
    filename = f"results_h_{h:0.4f}_{phasefield_model}"
    mesh_generation(0.1, h)
    # Convert GMSH (.msh) meshes to VTU meshes appropriate for OGS simulation.
    input_file = f"{out_dir}/" + meshname + ".msh"
    meshes = ot.meshes_from_gmsh(filename=input_file, log=False)
    for name, mesh in meshes.items():
        pv.save_meshio(f"{out_dir}/mesh_full_pf_{name}.vtu", mesh)
    %cd {out_dir}
    run(
        "identifySubdomains -f -m mesh_full_pf_domain.vtu -- mesh_full_pf_physical_group_*.vtu",
        shell=True,
        check=True,
    )
    %cd -

    # As a preprocessing step, define the initial phase-field (crack).
    pre_processing(h, a0)
    # change properties in prj file #For more information visit: https://ogstools.opengeosys.org/stable/reference/ogstools.ogs6py.html
    model = ot.Project(
        input_file=prj_name,
        output_file=f"{out_dir}/{prj_name}",
        MKL=True,
        args=f"-o {out_dir}",
    )
    model.replace_parameter_value(name="ls", value=ls)
    model.replace_text(phasefield_model, xpath="./processes/process/phasefield_model")
    model.replace_text(filename, xpath="./time_loop/output/prefix")
    model.write_input()
    # run simulation with ogs
    t0 = time.time()
    print(">>> OGS started execution ... <<<")
    run(
        f"ogs {out_dir}/{prj_name} -o {out_dir} > {out_dir}/log.txt",
        shell=True,
        check=True,
    )
    tf = time.time()
    print(">>> OGS terminated execution  <<< Elapsed time: ", round(tf - t0, 2), " s.")

def Analytical_solution(phasefield_model, h):
    v = np.linspace(1e-10, 0.3, 31)
    pres = np.linspace(0, 1.0, 31)
    length = np.linspace(0, 1.0, 31)

    ls = 2 * h

    # Effective Gc and a for AT1/A2
    if phasefield_model == "AT1":
        Gc_ref = Gc * (3 * h / 8.0 / ls + 1.0)
        a_eff = a0 * (1 + pi * ls / (4.0 * a0 * (3 * h / 8.0 / ls + 1.0)))
    elif phasefield_model == "AT2":
        Gc_ref = Gc * (h / (2.0 * ls) + 1.0)
        a_eff = a0 * (1 + pi * ls / (4.0 * a0 * (h / (2.0 * ls) + 1.0)))

    Eprime = E / (1 - nu**2)

    V_c = (4 * pi * Gc_ref * a_eff**3 / Eprime) ** 0.5
    P_c = (Eprime * Gc_ref / (pi * a_eff)) ** 0.5

    for i in range(len(v)):
        if v[i] < V_c:
            pres[i] = Eprime * v[i] / (2 * pi * a_eff**2) / P_c
            length[i] = a_eff
        else:
            pres[i] = (2 * Eprime * Gc_ref**2 / (pi * v[i])) ** 0.333333 / P_c
            length[i] = (Eprime * v[i] ** 2 / (4 * pi * Gc_ref)) ** 0.333333

    return v, pres, length, Gc_ref, P_c


fluidVolume_analytical = Analytical_solution(phasefield_model, h)[0]
pressure_analytical = Analytical_solution(phasefield_model, h)[1]
length_analytical = Analytical_solution(phasefield_model, h)[2]
Gc_ref = Analytical_solution(phasefield_model, h)[3]
P_c = Analytical_solution(phasefield_model, h)[4]

plt.subplots(figsize=(12, 4))
plt.subplot(1, 2, 1)
fluid_volume = []
surface_energy = []
# Open the file for reading
with (out_dir / "log.txt").open() as fd:
    # Iterate over the lines
    for _i, line in enumerate(fd):
        match_surface_energy = re.search(
            r"""Surface energy: (\d+\.\d*) Pressure work: (\d+\.\d*)  at time: (\d+\.\d*)""",
            line,
        )
        if match_surface_energy:
            surface_energy.append(float(match_surface_energy.group(1)))
            fluid_volume.append(float(match_surface_energy.group(3)))

plt.grid(linestyle="dashed")
plt.xlabel("Volume")
plt.ylabel("Crack length")
plt.plot(fluidVolume_analytical, length_analytical, "black", label="Closed form")
plt.plot(
    fluid_volume,
    np.array(surface_energy[:]) / Gc_ref / 2,
    "bo",
    fillstyle="none",
    label="Phase-field",
)
legend = plt.legend(loc="lower right")

plt.subplot(1, 2, 2)
fluid_volume = []
pressure = []
# Open the file for reading
with (out_dir / "log.txt").open() as fd:
    # Iterate over the lines
    for _i, line in enumerate(fd):
        match_pressure = re.search(
            r"""Pressure: (\d+\.\d*) at time: (\d+\.\d*)""", line
        )
        if match_pressure:
            fluid_volume.append(float(match_pressure.group(2)))
            pressure.append(float(match_pressure.group(1)))


plt.xlabel("Volume")
plt.ylabel("$p_f/p_c$")
plt.plot(fluidVolume_analytical, pressure_analytical, "black", label="Closed form")
plt.plot(
    np.array(fluid_volume[:]),
    np.array(pressure[:]) / P_c,
    "g<",
    fillstyle="none",
    label="Phase-field",
)
plt.grid(linestyle="dashed")
legend = plt.legend(loc="upper right")
plt.show()

plt.subplots(figsize=(12, 4))
plt.subplot(1, 2, 1)
plt.grid(linestyle="dashed")
plt.plot(
    fluidVolume_analytical[1:],
    (
        abs(length_analytical[1:] - np.array(surface_energy[:]) / (2 * Gc_ref))
        / (np.array(surface_energy[:]) / (2 * Gc_ref))
    )
    * 100,
    "-.ob",
    fillstyle="none",
    label="Closed form",
)
plt.ylim([0.0, 10])
plt.xlabel("Volume")
plt.ylabel(
    r"$\frac{|a_\mathrm{num}-{a}_\mathrm{ana}|}{{a}_\mathrm{num}}\times 100\%$",
    fontsize=14,
)

plt.subplot(1, 2, 2)
plt.grid(linestyle="dashed")
plt.plot(
    fluidVolume_analytical[1:],
    (abs(pressure_analytical[1:] - np.array(pressure[:]) / P_c) / pressure[:]) * 100,
    "-.<g",
    fillstyle="none",
    label="Closed form",
)
plt.ylim([0.0, 10])
plt.xlabel("Volume")
plt.ylabel(
    r"$\frac{|p_\mathrm{num}-{p}_\mathrm{ana}|}{{p}_\mathrm{num}}\times 100\%$",
    fontsize=14,
)
plt.show()

filename = f"results_h_{h:0.4f}_{phasefield_model}"
reader = pv.get_reader(f"{out_dir}/" + filename + ".pvd")

plotter = pv.Plotter(shape=(1, 2), border=False)
plotter.open_gif(f"{out_dir}/Kregime_Propagating.gif")
for time_value in reader.time_values:
    reader.set_active_time_value(time_value)
    mesh = reader.read()[0]
    sargs = {
        "title": "Phase field",
        "title_font_size": 16,
        "label_font_size": 12,
        "n_labels": 5,
        "position_x": 0.25,
        "position_y": 0.15,
        "fmt": "%.1f",
        "width": 0.5,
    }
    p = pv.Plotter(shape=(1, 2), border=False)
    clim = [0, 1.0]
    points = mesh.point_data["phasefield"].shape[0]
    xs = mesh.points[:, 0]
    ys = mesh.points[:, 1]
    pf = mesh.point_data["phasefield"]
    plotter.clear()
    plotter.add_mesh(
        mesh,
        scalars=pf,
        show_scalar_bar=True,
        colormap="coolwarm",
        clim=clim,
        scalar_bar_args=sargs,
        lighting=False,
    )
    plotter.add_text(f"Time: {time_value:0.02f}", color="black")

    plotter.view_xy()
    plotter.write_frame()

plotter.close()

mesh = reader.read()[0]

points = mesh.point_data["phasefield"].shape[0]
xs = mesh.points[:, 0]
ys = mesh.points[:, 1]
pf = mesh.point_data["phasefield"]

clim = [0, 1.0]
sargs = {
    "title": "Phase field",
    "title_font_size": 16,
    "label_font_size": 12,
    "n_labels": 5,
    "position_x": 0.25,
    "position_y": 0.0,
    "fmt": "%.1f",
    "width": 0.5,
}
plotter = pv.Plotter(shape=(1, 2), border=False)
plotter.add_mesh(
    mesh,
    scalars=pf,
    show_edges=False,
    show_scalar_bar=True,
    colormap="coolwarm",
    clim=clim,
    scalar_bar_args=sargs,
    lighting=False,
)

plotter.view_xy()
plotter.camera.zoom(1.5)
plotter.window_size = [1000, 500]
plotter.show()