# -*- coding: utf-8 -*-
"""Bigoni doubly clamped elastica: stable antisymmetric second mode."""

import sys

import numpy as np

sys.path.insert(0, "./server/")

from pytorsocae import TorsoCAESession


LENGTH = 0.490
WIDTH = 0.025
THICKNESS = 0.0015
YOUNGS_MODULUS = 3.0e9
POISSON_RATIO = 0.35
NORMALIZED_END_SHORTENING = 1.2
END_SHORTENING = NORMALIZED_END_SHORTENING * LENGTH
IMPERFECTION = 0.002 * LENGTH
MESH_DIVISIONS = (160, 4, 2)


def second_mode_seed(u, v, w):
    """Apply a small antisymmetric imperfection to select the second branch."""
    phase = 2.0 * np.pi * u
    z_seed = IMPERFECTION * np.sin(phase) * (1.0 - np.cos(phase))
    return (
        LENGTH * u,
        WIDTH * (v - 0.5),
        THICKNESS * (w - 0.5) + z_seed,
    )


def build_session():
    session = TorsoCAESession()
    session.build_structured_mesh(
        mesh_id="mesh_0",
        name="Bigoni Antisymmetric Elastica",
        body_id="body_0",
        coordinate_system="parametric",
        bounds=((0.0, 1.0), (0.0, 1.0), (0.0, 1.0)),
        divisions=MESH_DIVISIONS,
        surface_tags={
            "xmin": 1,
            "xmax": 2,
            "ymin": 3,
            "ymax": 4,
            "zmin": 5,
            "zmax": 6,
        },
        mapping=second_mode_seed,
    )
    session.solid("Bigoni Antisymmetric Elastica").material(
        E=YOUNGS_MODULUS,
        nu=POISSON_RATIO,
    )
    session.surface(1, scope_id="mesh_0").bc("fixed")
    session.surface(2, scope_id="mesh_0").bc(
        "fixed",
        values=[f"-{END_SHORTENING:.17g}*t", 0.0, 0.0],
    )
    session.set_physics(
        "structural",
        submodel="geometric_nonlinear",
        backend="dolfinx",
    )
    session.set_solver_options(
        algo="direct",
        precond="none",
        tol=1.0e-10,
        max_iter=500,
        num_steps=96,
        max_inner_iter=35,
        inner_tol=1.0e-8,
        max_load_bisections=10,
        n_cores=2,
        device="cpu",
    )
    return session


def run_benchmark():
    return build_session().compute(mesh_ids=["mesh_0"])


if __name__ == "__main__":
    result = run_benchmark()
    print(
        "Bigoni antisymmetric elastica:",
        f"max_displacement={result['max_displacement']:.6g} m",
        f"max_von_mises={result['max_von_mises']:.6g} Pa",
    )