# %%
# 2D lid-driven cavity flow with PFEM
# ===================================
# This example simulates a 2D lid-driven cavity flow using the Particle Finite Element Method (PFEM).
#
# %%
# Keywords
# --------
# PFEM, lid driven cavity, 2D, fluid
#
# %%
# Description
# -----------
# The test demonstrates:
# - How to perform a PFEM simulation in a fixed domain (without free surfaces)
import os, sys, shutil
import numpy as np

# %%
# Gmsh pfem branch
# ----------------
# This test requires the alphashapes branch of GMSH to use the Particle Finite Element (PFEM) module.
# PFEM is a lagrangian method where the mesh is regenerated at each time step, making it suitable for problems with large deformations and free surfaces.
sys.path.insert(0, os.environ["GMSH_PFEM_DIR"])

import gmsh
from migflow import fluid, pfem

# %%
# Output Directory
# ----------------
# Create a clean output directory for simulation results.
outputdir = "output_2d_lid_driven_cavity_pfem" if len(sys.argv) < 2 else sys.argv[1]
shutil.rmtree(outputdir, ignore_errors=True)
os.makedirs(outputdir)


# %%
# Geometrical parameters and mesh generation
# ------------------------------------------
l_box = 1.0
h_box = 1.0

# %%
# Mesh parameters
# ---------------
geo_mesh_size = l_box / 10
mesh_size = l_box / 40
alpha = 10.0

# %%
# Initial fluid mesh
# ------------------
gmsh.initialize()

gmsh.model.add("ModelInit")
gmsh.model.occ.addRectangle(0, 0, 0, l_box, h_box)
gmsh.model.occ.synchronize()
gmsh.model.mesh.setSizeCallback(lambda *args: mesh_size)
gmsh.model.mesh.generate(2)
nodeTags, coords, _ = gmsh.model.mesh.getNodes()

# %%
# Solid domain
# ------------
gmsh.model.add("ModelGeo")
gmsh.model.occ.addRectangle(0, 0, 0, l_box, h_box)
gmsh.model.occ.synchronize()
gmsh.model.mesh.setSizeCallback(lambda *args: mesh_size)
gmsh.model.mesh.generate(2)
nodeTags, coords, _ = gmsh.model.mesh.getNodes()
gmsh.model.addPhysicalGroup(1, [1], -1, "bottom")
gmsh.model.addPhysicalGroup(1, [2], -1, "right")
gmsh.model.addPhysicalGroup(1, [3], -1, "top")
gmsh.model.addPhysicalGroup(1, [4], -1, "left")
geoEntities = gmsh.model.getEntities(1)
# %%
# PFEM mesh and size fields
# -------------------------
gmsh.model.add("ModelFluid")
alphaDomainTag = gmsh.model.addDiscreteEntity(2, -1, [])
for dim, tag in geoEntities:
    gmsh.model.addDiscreteEntity(dim, tag, [])
alphaBoundaryTag = gmsh.model.addDiscreteEntity(1, -1, [])
gmsh.model.mesh.addNodes(2, alphaDomainTag, nodeTags, coords)
gmsh.model.addPhysicalGroup(1, [1], -1, "bottom")
gmsh.model.addPhysicalGroup(1, [2], -1, "right")
gmsh.model.addPhysicalGroup(1, [3], -1, "top")
gmsh.model.addPhysicalGroup(1, [4], -1, "left")
gmsh.model.addPhysicalGroup(1, [alphaBoundaryTag], -1, "freeSurface")
gmsh.model.addPhysicalGroup(2, [alphaDomainTag], -1, "domain")

# %%
# Size fields
# -----------
sizeFieldConstant = gmsh.model.mesh.field.add("Box")
gmsh.model.mesh.field.setNumber(sizeFieldConstant, "VIn", mesh_size)
gmsh.model.mesh.field.setNumber(sizeFieldConstant, "VOut", mesh_size)
gmsh.model.mesh.field.setNumber(sizeFieldConstant, "XMin", 0.0)
gmsh.model.mesh.field.setNumber(sizeFieldConstant, "XMax", l_box)
gmsh.model.mesh.field.setNumber(sizeFieldConstant, "YMin", 0.0)
gmsh.model.mesh.field.setNumber(sizeFieldConstant, "YMax", h_box)
gmsh.model.mesh.field.setNumber(sizeFieldConstant, "Thickness", 0.001)

# %%
# Physical Parameters
# -------------------
g = np.array([0.0, 0.0])
rho = 1000
mu = 1e-3
Re = 1000
v_top = Re * mu / (rho * l_box)
# %%
# Time parameters
# ---------------
cfl = 0.5
U = v_top
U_init = v_top
dt = mesh_size / U * cfl
t = 0
tEnd = 25000.0

# %%
# Fluid problem
# -------------
f = fluid.FluidProblem2(g, mu, rho, advection=False)
f.set_wall_boundary("bottom")
f.set_wall_boundary("right")
f.set_wall_boundary("top", velocity=[v_top, 0])
f.set_wall_boundary("left")
f.set_strong_boundary("bottom", velocity=[0, 0])
f.set_strong_boundary("right", velocity=[0, 0])
f.set_strong_boundary("top", velocity=[v_top, 0])
f.set_strong_boundary("left", velocity=[0, 0])

# %%
# Simulation Loop
# ---------------
# Time integration of coupled fluid–particle motion.
i = 0
outf = 25
gmsh.option.setNumber("General.Verbosity", 0)

while t < tEnd:
    print(f"{i:4d}, {t:.6g}/{tEnd:.6g}, {dt:.6g}")
    # Update PFEM mesh
    nodetag, oelemtag, oparamcoord, _ = gmsh.model.mesh.computeAlphaShape(
        2,
        alphaDomainTag,
        alphaBoundaryTag,
        "ModelGeo",
        alpha,
        sizeFieldConstant,
        sizeFieldConstant,
        boundaryTolerance=0.01 * mesh_size,
        usePreviousMesh=False,
    )
    oparamcoord = oparamcoord.reshape((-1, 3))
    gmsh.write(outputdir + "/lastMesh.msh")
    ordered_node_tags = pfem.prepareMeshForMigflow(
        i, alphaDomainTag, f, nodetag, oelemtag, oparamcoord
    )

    if i % outf == 0:  # or t > 0.3:
        f.write_mig(outputdir, t)

    # nodes velocity, be aware that its dimension is (n_nodes, 3) not (n_nodes, 2)
    u_old = np.zeros_like(f.coordinates())
    u_old[:, :2] = f.velocity()

    f.implicit_euler(dt)

    dx = np.zeros((f.coordinates().shape[0], 3))
    u = np.zeros_like(dx)
    u[:, :2] = f.velocity()
    if i == 0:
        u_old = u
    dx = u * dt + 0.5 * (u - u_old) / dt * dt**2

    # Do not move the top boundary nodes
    top_nodes = np.unique(f.mesh_boundaries()["top"].flatten())
    dx[top_nodes, :] = 0.0

    # advect nodes and project if needed
    gmsh.model.mesh.advectMeshNodes(
        2,
        alphaDomainTag,
        alphaBoundaryTag,
        "ModelGeo",
        ordered_node_tags,
        dx.flatten(),
        0.01 * mesh_size,
    )

    # set new coordinates
    _, newCoords, _ = gmsh.model.mesh.getNodes(2, alphaDomainTag)
    f.set_coordinates(newCoords)

    t += dt
    i += 1

# %%
# Plot
# ----
# .. code-block:: shell
#
#  python3 -m migflow.plot.migplot output_2d_lid_driven_cavity_pfem --actors fluid --fluid-field velocity --fluid-vmin 0 --fluid-vmax 0.001 --show-edges 1
# %%
# Artifacts
# ---------
# - output_2d_lid_driven_cavity_pfem/animation.mp4