"""Ishihara et al. (2023), P-15 3D bimorph, Table 6, 3 mm width.

DOI 10.32604/cmes.2023.024614. The CAD, Hex27 mesh, material, electrodes,
supports, solve, and tip measurement all use TorsoCAE session operations.
The paper used 960 Hex20 cells; the same 40 x 6 x 4 divisions are built here
with complete quadratic Hex27 cells. All dimensions and material data are SI.
"""
import numpy as np
from pytorsocae import TorsoCAESession


LENGTH = 0.020
WIDTH = 0.003
LAYER_THICKNESS = 0.00025
DRIVE_VOLTAGE = 1.0
REFERENCE_TIP_UM = 0.79706  # Table 6, P-15, w = 3 mm, strong coupling
MESH_ID = "ishihara_p15_mesh"
BODY_ID = "ishihara_p15_bimorph"

session = TorsoCAESession()
full = session.csg_make_shape(
    "box", {"dx": LENGTH, "dy": WIDTH, "dz": 2 * LAYER_THICKNESS,
            "cx": 0, "cy": 0, "cz": -LAYER_THICKNESS},
    name="Bimorph stock", body_id="bimorph_stock",
)
upper_tool = session.csg_make_shape(
    "box", {"dx": LENGTH, "dy": WIDTH, "dz": LAYER_THICKNESS,
            "cx": 0, "cy": 0, "cz": 0},
    name="Electrode partition tool", body_id="electrode_tool",
)
halves = session.csg_boolean(
    "partition", full, upper_tool, name="Two piezoelectric layers",
    body_id="bimorph_halves",
)
middle_tool = session.csg_make_shape(
    "box", {"dx": LENGTH, "dy": WIDTH, "dz": LAYER_THICKNESS,
            "cx": 0, "cy": 0, "cz": -LAYER_THICKNESS / 2},
    name="Mesh partition tool", body_id="mesh_partition_tool",
)
body = session.csg_boolean(
    "partition", halves, middle_tool, name="P-15 piezoelectric bimorph",
    body_id=BODY_ID,
)
for intermediate in (middle_tool, halves, upper_tool, full):
    session.csg_delete(intermediate)
session.csg_select(body)
mesh = session.mesh(
    algo_id="structured_hex", size_max=0.0005, order=2, mesh_id=MESH_ID,
)
assert mesh["element_types"].get("Hex27") == 960, "Expected 40 x 6 x 4 Hex27 cells"
assert mesh["topology"]["n_volumes"] == 4, "The CAD partitions were not retained"

session.set_physics("electromagnetics", "piezoelectric", backend="dolfinx")
boundaries = session.mesh_boundaries(MESH_ID)


def on_plane(boundary, axis, value):
    box = boundary["bbox"]
    return (abs(box[f"{axis}min"] - value) < 1e-11
            and abs(box[f"{axis}max"] - value) < 1e-11)


center = [face["tag"] for face in boundaries if on_plane(face, "z", 0)]
outer = [face["tag"] for face in boundaries
         if on_plane(face, "z", -LAYER_THICKNESS)
         or on_plane(face, "z", LAYER_THICKNESS)]
root = [face["tag"] for face in boundaries if on_plane(face, "x", 0)]
assert len(center) == 1 and len(outer) == 2 and len(root) == 4

# Mesh scope is deliberate: the paper's two line supports select mesh nodes,
# while the center electrode is an internal mesh facet. Re-meshing requires
# rebuilding these selections by replaying this journal.
for tag in center:
    session.surface(tag, scope_id=MESH_ID).bc(
        "potential", [DRIVE_VOLTAGE], physics="electromagnetics",
    )
for tag in outer:
    session.surface(tag, scope_id=MESH_ID).bc("ground", physics="electromagnetics")
for tag in root:
    session.surface(tag, scope_id=MESH_ID).bc(
        "fixed", components=["x"], physics="structural",
    )

# Figure 5, Case A: line AD holds y, line CD holds z at the left end.
points = np.asarray(session.mesh_data(MESH_ID)["x"], dtype=float)
on_root = np.isclose(points[:, 0], 0, rtol=0, atol=1e-12)
line_ad = np.flatnonzero(
    on_root & np.isclose(points[:, 1], WIDTH, rtol=0, atol=1e-12)
).tolist()
line_cd = np.flatnonzero(
    on_root & np.isclose(points[:, 2], -LAYER_THICKNESS, rtol=0, atol=1e-12)
).tolist()
assert line_ad and line_cd, "The root edge nodes were not found"
session.surface(-1, scope_id=MESH_ID).bc(
    "fixed", components=["y"], point_node_indices=line_ad,
    physics="structural",
)
session.surface(-2, scope_id=MESH_ID).bc(
    "fixed", components=["z"], point_node_indices=line_cd,
    physics="structural",
)

# P-15 at constant E, Table 1. Voigt order is xx, yy, zz, yz, xz, xy.
c11, c12, c13, c33, c44, c66 = [
    value * 1e10 for value in (13.7, 8.81, 8.83, 12.0, 1.71, 2.44)
]
stiffness = [
    [c11, c12, c13, 0, 0, 0],
    [c12, c11, c13, 0, 0, 0],
    [c13, c13, c33, 0, 0, 0],
    [0, 0, 0, c44, 0, 0],
    [0, 0, 0, 0, c44, 0],
    [0, 0, 0, 0, 0, c66],
]
coupling = [
    [0, 0, 0, 0, 14.4, 0],
    [0, 0, 0, 14.4, 0, 0],
    [-10.8, -10.8, 8.7, 0, 0, 0],
]
vacuum_permittivity = 8.8541878128e-12
relative_permittivity = [
    [1.90e-8 / vacuum_permittivity, 0, 0],
    [0, 1.90e-8 / vacuum_permittivity, 0],
    [0, 0, 1.24e-8 / vacuum_permittivity],
]
for name in session.mesh_data(MESH_ID)["volume_names"].values():
    session.solid(name).material(
        physics="piezoelectric", elastic_stiffness=stiffness,
        piezoelectric_coupling=coupling,
        relative_permittivity=relative_permittivity,
    )

session.set_solver_options(algo="direct", device="cpu", n_cores=2, tol=1e-10)
assert session.check_supports()["status"] == "ok"
result = session.compute()
profile = session.report(
    "line_profile", "disp_z", result=result,
    line_start=[0, WIDTH / 2, 0], line_end=[LENGTH, WIDTH / 2, 0],
    n_points=41, output_dir="",
)
tip_um = float(profile["values"][-1]) * 1e6
error_percent = 100 * (tip_um - REFERENCE_TIP_UM) / REFERENCE_TIP_UM
print(f"P-15 tip center: Torso {tip_um:.8f} µm; paper {REFERENCE_TIP_UM:.5f} µm; "
      f"error {error_percent:+.4f}%")
if abs(error_percent) > 5:
    raise AssertionError("Ishihara Table 6 tip displacement differs by more than 5%")