#!/usr/bin/env bash
# meshDictCreator — interactive cfMesh meshDict generator with box/sphere/(semi)cylinder enclosures
# Requires: OpenFOAM tools (surfaceCheck, surfaceGenerateBoundingBox, surfaceFeatureEdges,
#           surfaceTransformPoints, blockMesh, foamToSurface, patchview)

set -euo pipefail

red()   { printf "\033[31m%s\033[0m\n" "$*"; }
green() { printf "\033[32m%s\033[0m\n" "$*"; }
cyan()  { printf "\033[36m%s\033[0m\n" "$*"; }
bold()  { printf "\033[1m%s\033[0m\n"   "$*"; }

ask() {
  local p="$1"; local d="${2-}"
  if [[ -n "$d" ]]; then read -r -p "$p [$d]: " ans || true; echo "${ans:-$d}";
  else read -r -p "$p: " ans || true; echo "$ans"; fi
}

# For meshDict values: allow 'not'
ensure_number() { local v="$1"; [[ "$v" == "not" ]] && return 0; [[ "$v" =~ ^-?[0-9]*\.?[0-9]+$ ]] || { red "Expected a number: '$v'"; exit 1; }; }
# For mandatory numbers (bbox, angle): do NOT allow 'not'
ensure_number_strict() { local v="$1"; [[ "$v" =~ ^-?[0-9]*\.?[0-9]+$ ]] || { red "Expected a number: '$v'"; exit 1; }; }

# --- Patchview cache cleanup on exit (Windows/WSL) ---
cleanup_patchview_cache() {
  if command -v powershell.exe >/dev/null 2>&1; then
    powershell.exe -NoLogo -NoProfile -Command "Remove-Item -Path 'C:\pvtmp\patchview\*' -Force -Recurse -ErrorAction SilentlyContinue" || true
  fi
  if [[ -d /mnt/c/pvtmp/patchview ]]; then
    rm -rf /mnt/c/pvtmp/patchview/* 2>/dev/null || true
  fi
}
trap cleanup_patchview_cache EXIT

# Globals
ENC_DIR="/home/deadbacco/OpenFOAM/deadbacco-v2306/platforms/linux64GccDPInt32Opt/bin/enclosureGeometries"

# --- Patch preview helpers (ParaView) ---
list_patches() {
  echo "Patches (from combined STL):"
  printf ' - %s\n' "${patches_unique[@]}"
  echo "Type 'parapatch' to preview in ParaView (runs: patchview ${COMB})."
}


run_patchview() {
  local stl="$1"
  # try PATH first
  if command -v patchview >/dev/null 2>&1; then
    patchview "$stl" & disown || true
    return 0
  fi
  # try FOAM_USER_APPBIN
  if [[ -n "${FOAM_USER_APPBIN-}" && -x "$FOAM_USER_APPBIN/patchview" ]]; then
    "$FOAM_USER_APPBIN/patchview" "$stl" & disown || true
    return 0
  fi
  # try common OpenFOAM user platforms
  local cand
  cand="$(compgen -G "$HOME/OpenFOAM/${USER}-*/platforms/*/bin/patchview" | head -n1 || true)"
  if [[ -n "$cand" && -x "$cand" ]]; then
    "$cand" "$stl" & disown || true
    return 0
  fi
  red "patchview not found. Put it in PATH or \$FOAM_USER_APPBIN, then retry 'parapatch'."
  return 1
}


maybe_run_parapatch() {
  local v="$1"
  if [[ "$v" == "parapatch" ]]; then
    if [[ -f "$COMB" ]]; then
      cyan "Opening patch preview: patchview '$COMB'"
      # Interaktif login-shell, hataları göster, ÖN PLANDA çalışsın
      bash -lic "patchview '$COMB'"
    else
      red "Combined STL ($COMB) not found yet."
    fi
    return 0
  fi
  return 1
}





# ---- Cylinder blockMesh writers (based on your templates) ----

# Full cylinder (8 wedge blocks). Axis = Z. Center at (CX,CY,CZ). Radius=R. Height=H.
# Cell counts: nWall (radial per wedge), nCirc (around), nAx (axial).
write_cylinder_bmd_full() {
  local dictPath="$1"; shift
  local R="$1" H="$2" CX="$3" CY="$4" CZ="$5" nWall="$6" nCirc="$7" nAx="$8"

  local z0; local z1
  z0=$(awk "BEGIN{printf \"%.12f\", $CZ - $H/2.0}")
  z1=$(awk "BEGIN{printf \"%.12f\", $CZ + $H/2.0}")

  # unit circle points (9 around)
  local x1=1      y1=0
  local x2=0.7071 y2=0.7071
  local x3=0      y3=1
  local x4=-0.7071 y4=0.7071
  local x5=-1     y5=0
  local x6=-0.7071 y6=-0.7071
  local x7=0      y7=-1
  local x8=0.7071 y8=-0.7071

  # helper scale/shift
  ss(){ awk "BEGIN{printf \"%.12f\", $1*$2 + $3}"; }

  cat >"$dictPath" <<EOF
/*--------------------------------*- C++ -*----------------------------------*\\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2306                                  |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
convertToMeters 1.0;

vertices
(
    ($(ss $R 0 $CX) $(ss $R 0 $CY) $z0)   // 0 center bottom
    ($(ss $R $x1 $CX) $(ss $R $y1 $CY) $z0)   // 1
    ($(ss $R $x2 $CX) $(ss $R $y2 $CY) $z0)   // 2
    ($(ss $R $x3 $CX) $(ss $R $y3 $CY) $z0)   // 3
    ($(ss $R $x4 $CX) $(ss $R $y4 $CY) $z0)   // 4
    ($(ss $R $x5 $CX) $(ss $R $y5 $CY) $z0)   // 5
    ($(ss $R $x6 $CX) $(ss $R $y6 $CY) $z0)   // 6
    ($(ss $R $x7 $CX) $(ss $R $y7 $CY) $z0)   // 7
    ($(ss $R $x8 $CX) $(ss $R $y8 $CY) $z0)   // 8

    ($(ss $R 0 $CX) $(ss $R 0 $CY) $z1)   // 9 center top
    ($(ss $R $x1 $CX) $(ss $R $y1 $CY) $z1)   //10
    ($(ss $R $x2 $CX) $(ss $R $y2 $CY) $z1)   //11
    ($(ss $R $x3 $CX) $(ss $R $y3 $CY) $z1)   //12
    ($(ss $R $x4 $CX) $(ss $R $y4 $CY) $z1)   //13
    ($(ss $R $x5 $CX) $(ss $R $y5 $CY) $z1)   //14
    ($(ss $R $x6 $CX) $(ss $R $y6 $CY) $z1)   //15
    ($(ss $R $x7 $CX) $(ss $R $y7 $CY) $z1)   //16
    ($(ss $R $x8 $CX) $(ss $R $y8 $CY) $z1)   //17
);

blocks
(
    hex (0 1 2 0 9 10 11 9) ($nWall $nCirc $nAx) simpleGrading (1 1 1)
    hex (0 2 3 0 9 11 12 9) ($nWall $nCirc $nAx) simpleGrading (1 1 1)
    hex (0 3 4 0 9 12 13 9) ($nWall $nCirc $nAx) simpleGrading (1 1 1)
    hex (0 4 5 0 9 13 14 9) ($nWall $nCirc $nAx) simpleGrading (1 1 1)
    hex (0 5 6 0 9 14 15 9) ($nWall $nCirc $nAx) simpleGrading (1 1 1)
    hex (0 6 7 0 9 15 16 9) ($nWall $nCirc $nAx) simpleGrading (1 1 1)
    hex (0 7 8 0 9 16 17 9) ($nWall $nCirc $nAx) simpleGrading (1 1 1)
    hex (0 8 1 0 9 17 10 9) ($nWall $nCirc $nAx) simpleGrading (1 1 1)
);

edges
(
    arc 1 2 ($(ss $R 0.9239 $CX) $(ss $R 0.3827 $CY) $z0)
    arc 2 3 ($(ss $R 0.3827 $CX) $(ss $R 0.9239 $CY) $z0)
    arc 3 4 ($(ss $R -0.3827 $CX) $(ss $R 0.9239 $CY) $z0)
    arc 4 5 ($(ss $R -0.9239 $CX) $(ss $R 0.3827 $CY) $z0)
    arc 5 6 ($(ss $R -0.9239 $CX) $(ss $R -0.3827 $CY) $z0)
    arc 6 7 ($(ss $R -0.3827 $CX) $(ss $R -0.9239 $CY) $z0)
    arc 7 8 ($(ss $R 0.3827 $CX) $(ss $R -0.9239 $CY) $z0)
    arc 8 1 ($(ss $R 0.9239 $CX) $(ss $R -0.3827 $CY) $z0)

    arc 10 11 ($(ss $R 0.9239 $CX) $(ss $R 0.3827 $CY) $z1)
    arc 11 12 ($(ss $R 0.3827 $CX) $(ss $R 0.9239 $CY) $z1)
    arc 12 13 ($(ss $R -0.3827 $CX) $(ss $R 0.9239 $CY) $z1)
    arc 13 14 ($(ss $R -0.9239 $CX) $(ss $R 0.3827 $CY) $z1)
    arc 14 15 ($(ss $R -0.9239 $CX) $(ss $R -0.3827 $CY) $z1)
    arc 15 16 ($(ss $R -0.3827 $CX) $(ss $R -0.9239 $CY) $z1)
    arc 16 17 ($(ss $R 0.3827 $CX) $(ss $R -0.9239 $CY) $z1)
    arc 17 10 ($(ss $R 0.9239 $CX) $(ss $R -0.3827 $CY) $z1)
);

boundary
(
    wall
    {
        type wall;
        faces
        (
            (1 2 11 10)
            (2 3 12 11)
            (3 4 13 12)
            (4 5 14 13)
            (5 6 15 14)
            (6 7 16 15)
            (7 8 17 16)
            (8 1 10 17)
        );
    }

    top
    {
        type patch;
        faces
        (
            (9 10 11 9)
            (9 11 12 9)
            (9 12 13 9)
            (9 13 14 9)
            (9 14 15 9)
            (9 15 16 9)
            (9 16 17 9)
            (9 17 10 9)
        );
    }

    bottom
    {
        type patch;
        faces
        (
            (0 1 2 0)
            (0 2 3 0)
            (0 3 4 0)
            (0 4 5 0)
            (0 5 6 0)
            (0 6 7 0)
            (0 7 8 0)
            (0 8 1 0)
        );
    }
);

mergePatchPairs();
EOF
}

# Semi-cylinder (half). Default flat face at y=0 plane (xz-plane). Axis = Z.
# If symAxis=x, we rotate exported STL by -90° about Z to make flat on x=0 plane.
write_semicyl_bmd_full() {
  local dictPath="$1"; shift
  local R="$1" H="$2" CX="$3" CY="$4" CZ="$5"

  local z0; local z1
  z0=$(awk "BEGIN{printf \"%.12f\", $CZ - $H/2.0}")
  z1=$(awk "BEGIN{printf \"%.12f\", $CZ + $H/2.0}")

  ss(){ awk "BEGIN{printf \"%.12f\", $1*$2 + $3}"; }

  cat >"$dictPath" <<EOF
/* system/blockMeshDict (semi-cylinder) */
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
convertToMeters 1;

vertices
(
    // z = $z0
    ($(ss $R -1 $CX)     $(ss $R 0 $CY)   $z0)   // 0
    ($(ss $R  0 $CX)     $(ss $R 0 $CY)   $z0)   // 1
    ($(ss $R  1 $CX)     $(ss $R 0 $CY)   $z0)   // 2
    ($(ss $R  0 $CX)     $(ss $R 1 $CY)   $z0)   // 3
    ($(ss $R -0.707106 $CX) $(ss $R 0.707106 $CY) $z0)   // 4
    ($(ss $R  0.707106 $CX) $(ss $R 0.707106 $CY) $z0)   // 5

    // z = $z1
    ($(ss $R -1 $CX)     $(ss $R 0 $CY)   $z1)   // 6
    ($(ss $R  0 $CX)     $(ss $R 0 $CY)   $z1)   // 7
    ($(ss $R  1 $CX)     $(ss $R 0 $CY)   $z1)   // 8
    ($(ss $R  0 $CX)     $(ss $R 1 $CY)   $z1)   // 9
    ($(ss $R -0.707106 $CX) $(ss $R 0.707106 $CY) $z1)   //10
    ($(ss $R  0.707106 $CX) $(ss $R 0.707106 $CY) $z1)   //11
);

blocks
(
    hex (0 1 3 4 6 7 9 10) (20 20 40) simpleGrading (1 1 1)
    hex (1 2 5 3 7 8 11 9) (20 20 40) simpleGrading (1 1 1)
);

edges
(
    // z=$z0 outer arcs
    arc 0 4 ($(ss $R -0.923880 $CX) $(ss $R 0.382683 $CY) $z0)
    arc 4 3 ($(ss $R -0.382683 $CX) $(ss $R 0.923880 $CY) $z0)
    arc 2 5 ($(ss $R 0.923880 $CX) $(ss $R 0.382683 $CY) $z0)
    arc 5 3 ($(ss $R 0.382683 $CX) $(ss $R 0.923880 $CY) $z0)
    // z=$z1 outer arcs
    arc 6 10 ($(ss $R -0.923880 $CX) $(ss $R 0.382683 $CY) $z1)
    arc 10 9 ($(ss $R -0.382683 $CX) $(ss $R 0.923880 $CY) $z1)
    arc 8 11 ($(ss $R 0.923880 $CX) $(ss $R 0.382683 $CY) $z1)
    arc 11 9 ($(ss $R 0.382683 $CX) $(ss $R 0.923880 $CY) $z1)
);

boundary
(
    curvedWall
    {
        type wall;
        faces
        (
            (0 4 10 6)
            (4 3 9 10)
            (2 5 11 8)
            (5 3 9 11)
        );
    }
    symm
    {
        type symmetryPlane;
        faces
        (
            (0 1 7 6)
            (1 2 8 7)
        );
    }
    bottom
    {
        type patch;
        faces
        (
            (0 1 3 4)
            (1 2 5 3)
        );
    }
    top
    {
        type patch;
        faces
        (
            (6 7 9 10)
            (7 8 11 9)
        );
    }
);

mergePatchPairs();
EOF
}

# 0) Pick STL (body)
shopt -s nullglob; stls=( *.stl *.STL ); shopt -u nullglob
(( ${#stls[@]} )) || { red "No .stl found."; exit 1; }
cand=()
for f in "${stls[@]}"; do
  case "$(basename "$f")" in
    sphere.stl|semisphere.stl|boundary.stl|*Box.stl) ;; # skip enclosures
    *) cand+=("$f");;
  esac
done
(( ${#cand[@]} )) || { red "No body STL candidate found."; exit 1; }
if (( ${#cand[@]} == 1 )); then STL="${cand[0]}"; else bold "STL candidates:"; printf ' - %s\n' "${cand[@]}"; STL="$(ask "Select body STL" "${cand[0]}")"; fi
[[ -f "$STL" ]] || { red "File not found: $STL"; exit 1; }
BASE="${STL%.*}"
COMB="${BASE}Box.stl"   # combined STL name derived from body

# 1) surfaceCheck
cyan "Running surfaceCheck..."
if ! surfaceCheck "$STL" >/tmp/meshDictCreator_surfaceCheck.log 2>&1; then
  red "surfaceCheck failed. See /tmp/meshDictCreator_surfaceCheck.log"; exit 1; fi
green "STL OK."

# 2) Symmetry + domain selection
use_sym="$(ask "Is there a symmetry plane? (y/n)" "n")"
sym_axis=""; sym_plane_axis=""
case "$use_sym" in
  [Yy]*) 
    sym_axis="$(ask "Body direction? (+x,-x,+y,-y,+z,-z)" "+x")"
    case "$sym_axis" in
      +x|-x) sym_plane_axis="x" ;;
      +y|-y) sym_plane_axis="y" ;;
      +z|-z) sym_plane_axis="z" ;;
      *) red "Invalid choice."; exit 1 ;;
    esac
    domain="$(ask "Choose enclosure: box / semisphere / semicylinder" "semisphere")"
    ;;
  *)
    domain="$(ask "Choose enclosure: box / sphere / cylinder" "box")"
    ;;
esac

# 3) Build enclosure -> produce $COMB
case "$domain" in
  box)
    bold "Bounding box limits (xMin xMax yMin yMax zMin zMax)."
    declare -A lim; declare -A fixed
    fix_plane(){ case "$1" in +x|-x) lim[xMin]=0; fixed[xMin]=1;; +y|-y) lim[yMin]=0; fixed[yMin]=1;; +z|-z) lim[zMin]=0; fixed[zMin]=1;; esac; }
    [[ -n "$sym_axis" ]] && fix_plane "$sym_axis"
    for k in xMin xMax yMin yMax zMin zMax; do
      if [[ -n "${fixed[$k]-}" ]]; then echo "$k = 0 (symmetry)"; else
        val="$(ask "$k value" "9")"; [[ "$val" == "not" ]] && { red "'not' not allowed for bounding box."; exit 1; }
        ensure_number_strict "$val"; lim[$k]="$val"; fi
    done
    cyan "Generating box enclosure STL..."
    surfaceGenerateBoundingBox "$STL" "$COMB" "${lim[xMin]}" "${lim[xMax]}" "${lim[yMin]}" "${lim[yMax]}" "${lim[zMin]}" "${lim[zMax]}"
    ;;

  sphere)
    # radius + sphere center offsets (sphere moves; body stays)
    R="$(ask "Sphere radius [m]" "5")"; ensure_number_strict "$R"
    dx="$(ask "Sphere center offset dx" "0")"; dy="$(ask "dy" "0")"; dz="$(ask "dz" "0")"
    ensure_number_strict "$dx"; ensure_number_strict "$dy"; ensure_number_strict "$dz"
    cp "$ENC_DIR/sphere.stl" ./sphere.stl
    cyan "Scaling sphere.stl..."
    surfaceTransformPoints -scale "($R $R $R)" sphere.stl sphere.tmp.stl
    mv -f sphere.tmp.stl sphere.stl
    cyan "Translating sphere.stl..."
    surfaceTransformPoints -translate "($dx $dy $dz)" sphere.stl sphere.tmp.stl
    mv -f sphere.tmp.stl sphere.stl
    cyan "Combining STL..."
    cat "$STL" sphere.stl > "$COMB"
    ;;

  semisphere)
    [[ -n "$sym_plane_axis" ]] || { red "Semi-sphere requires symmetry plane."; exit 1; }
    R="$(ask "Semi-sphere radius [m]" "5")"; ensure_number_strict "$R"
    cp "$ENC_DIR/semisphere.stl" ./semisphere.stl
    # Base semisphere: flat on z=0 (dome +z). Rotate with rollPitchYaw (deg) to align flat face to requested plane.
    case "$sym_plane_axis" in
      x) cyan "Orienting semi-sphere: flat -> x=0"; surfaceTransformPoints -rollPitchYaw "(0 90 0)" semisphere.stl semisphere.rot.stl ;;
      y) cyan "Orienting semi-sphere: flat -> y=0"; surfaceTransformPoints -rollPitchYaw "(-90 0 0)" semisphere.stl semisphere.rot.stl ;;
      z) cyan "Orienting semi-sphere: flat already at z=0"; cp semisphere.stl semisphere.rot.stl ;;
    esac
    mv -f semisphere.rot.stl semisphere.stl
    cyan "Scaling semi-sphere..."
    surfaceTransformPoints -scale "($R $R $R)" semisphere.stl semisphere.tmp.stl
    mv -f semisphere.tmp.stl semisphere.stl
    cyan "Combining STL..."
    cat "$STL" semisphere.stl > "$COMB"
    ;;

    cylinder)
    R="$(ask "Cylinder radius [m]" "5")"; ensure_number_strict "$R"
    H="$(ask "Cylinder height [m]" "10")"; ensure_number_strict "$H"
    CX="$(ask "Cylinder center X" "0")"; CY="$(ask "Y" "0")"; CZ="$(ask "Z" "0")"
    ensure_number_strict "$CX"; ensure_number_strict "$CY"; ensure_number_strict "$CZ"
    nWall="$(ask "cells radial (per wedge)" "10")"; nCirc="$(ask "cells circumferential" "10")"; nAx="$(ask "cells axial" "20")"
    ensure_number_strict "$nWall"; ensure_number_strict "$nCirc"; ensure_number_strict "$nAx"

    # NEW: ask flow axis
    axis="$(ask "Flow axis for cylinder (x/y/z)" "z")"

    mkdir -p system
    dictCyl="system/blockMeshDict.cyl"
    cyan "Writing blockMeshDict (full cylinder, default axis = z)..."
    write_cylinder_bmd_full "$dictCyl" "$R" "$H" "$CX" "$CY" "$CZ" "$nWall" "$nCirc" "$nAx"

    cyan "Running blockMesh..."
    blockMesh -dict "$dictCyl" >/tmp/meshDictCreator_blockMesh.log 2>&1

    cyan "Exporting boundary to STL (current dir)..."
    ( foamToSurface boundary.stl > /tmp/meshDictCreator_foamToSurface.log 2>&1 ) || true
    [[ -f ./boundary.stl ]] || { red "boundary.stl not produced."; exit 1; }

    # Rotate around cylinder center so its axis aligns with requested flow axis
    case "${axis,,}" in
      x)
        # z-axis -> x-axis: rotate -90° about Y, around (CX,CY,CZ)
        surfaceTransformPoints -translate "($(( -1 * CX )) $(( -1 * CY )) $(( -1 * CZ )))" boundary.stl boundary.tmp1.stl
        surfaceTransformPoints -rollPitchYaw "(0 -90 0)"        boundary.tmp1.stl boundary.tmp2.stl
        surfaceTransformPoints -translate  "($CX $CY $CZ)"      boundary.tmp2.stl boundary.stl
        rm -f boundary.tmp1.stl boundary.tmp2.stl
        ;;
      y)
        # z-axis -> y-axis: rotate +90° about X, around (CX,CY,CZ)
        surfaceTransformPoints -translate "($(( -1 * CX )) $(( -1 * CY )) $(( -1 * CZ )))" boundary.stl boundary.tmp1.stl
        surfaceTransformPoints -rollPitchYaw "(90 0 0)"         boundary.tmp1.stl boundary.tmp2.stl
        surfaceTransformPoints -translate  "($CX $CY $CZ)"      boundary.tmp2.stl boundary.stl
        rm -f boundary.tmp1.stl boundary.tmp2.stl
        ;;
      z|*)
        # already along z, no rotation
        :
        ;;
    esac

    cyan "Combining STL..."
    cat "$STL" boundary.stl > "$COMB"
    ;;


  semicylinder)
    [[ -n "$sym_plane_axis" ]] || { red "Semi-cylinder requires symmetry plane."; exit 1; }
    
    R="$(ask "Semi-cylinder radius [m]" "5")"; ensure_number_strict "$R"
    H="$(ask "Semi-cylinder height [m]" "10")"; ensure_number_strict "$H"
    CX="$(ask "Center X" "0")"; CY="$(ask "Y" "0")"; CZ="$(ask "Z" "0")"
    flow="$(ask "Flow direction? (+x,-x,+y,-y,+z,-z)" "+z")"
case "$flow" in +x|-x|+y|-y|+z|-z) ;; *) red "Invalid flow"; exit 1;; esac

    ensure_number_strict "$CX"; ensure_number_strict "$CY"; ensure_number_strict "$CZ"
    mkdir -p system
    dictCyl="system/blockMeshDict.semicyl"
    cyan "Writing blockMeshDict (semi-cylinder, flat at y=0)..."
    write_semicyl_bmd_full "$dictCyl" "$R" "$H" "$CX" "$CY" "$CZ"
    cyan "Running blockMesh..."
    blockMesh -dict "$dictCyl" >/tmp/meshDictCreator_blockMesh.log 2>&1
    cyan "Exporting boundary to STL (current dir)..."
    ( foamToSurface boundary.stl > /tmp/meshDictCreator_foamToSurface.log 2>&1 ) || true
    [[ -f ./boundary.stl ]] || { red "boundary.stl not produced."; exit 1; }

    # safe translate vector strings (negatif değerlerde --1 hatasını önlemek için)
neg_vec() {
  awk -v x="$1" -v y="$2" -v z="$3" 'BEGIN{printf "(%.12f %.12f %.12f)", -x, -y, -z}'
}
pos_vec() {
  awk -v x="$1" -v y="$2" -v z="$3" 'BEGIN{printf "(%.12f %.12f %.12f)",  x,  y,  z}'
}

center_neg="$(neg_vec "$CX" "$CY" "$CZ")"
center_pos="$(pos_vec "$CX" "$CY" "$CZ")"

# 3.1) EKSEN HİZALAMA: başlangıçta eksen Z. Z -> flow
# rollPitchYaw = (roll about X, pitch about Y, yaw about Z) [deg]
case "$flow" in
  +z) r1="(0 0 0)" ;;
  -z) r1="(180 0 0)" ;;          # Z -> -Z (X etrafında 180 de yeter)
  +x) r1="(0 -90 0)" ;;          # Z -> +X  (Y etrafında -90)
  -x) r1="(0 90 0)" ;;           # Z -> -X  (Y etrafında +90)
  +y) r1="(90 0 0)" ;;           # Z -> +Y  (X etrafında +90)
  -y) r1="(-90 0 0)" ;;          # Z -> -Y  (X etrafında -90)
esac

surfaceTransformPoints -translate "$center_neg" boundary.stl b_aln1.stl
surfaceTransformPoints -rollPitchYaw "$r1"      b_aln1.stl  b_aln2.stl
rm -f b_aln1.stl

# 3.2) DÜZ YÜZEY → SİMETRİ DÜZLEMİ
# Şablon yarım silindirde düz yüzey başlangıçta Y=0 (düzlemin normali +Y).
# Eksen akışa hizalandıktan sonra, düz yüzeyi istenen düzleme taşımak için
# YALNIZCA "akış ekseni etrafında" ikinci bir dönüş yap.
# Bu yüzden ikinci dönüşte sadece akış eksenine karşılık gelen açıyı kullan.
case "$flow" in
  +x|-x)  # akış ekseni X → ikinci dönüş X etrafında (roll)
    case "$sym_plane_axis" in
      x) r2="(0 0 0)" ;;          # düz yüzey x=0 zaten eksen düzleminde (çevirmeye gerek yok)
      y) r2="(0 0 0)" ;;          # başlangıçta Y=0 idi; eksen hizalaması Y’yi korur → değişme yok
      z) r2="(90 0 0)" ;;         # Y normalini Z’ye döndürmek için X etrafında +90
    esac
  ;;
  +y|-y)  # akış ekseni Y → ikinci dönüş Y etrafında (pitch)
    case "$sym_plane_axis" in
      y) r2="(0 0 0)" ;;
      x) r2="(0 -90 0)" ;;        # Y normalini X’e döndürmek için Z→X eşdeğer pitch -90
      z) r2="(0 90 0)" ;;         # Y normalini Z’ye döndürmek için pitch +90
    esac
  ;;
  +z|-z)  # akış ekseni Z → ikinci dönüş Z etrafında (yaw)
    case "$sym_plane_axis" in
      z) r2="(0 0 0)" ;;          
      x) r2="(0 0 -90)" ;;        # Y normalini X’e döndürmek için Z etrafında -90
      y) r2="(0 0 0)" ;;          # zaten Y=0 düzlemi
    esac
  ;;
esac

surfaceTransformPoints -rollPitchYaw "$r2"      b_aln2.stl  b_aln3.stl
surfaceTransformPoints -translate "$center_pos" b_aln3.stl  boundary.stl
rm -f b_aln2.stl b_aln3.stl


    # Align flat face to requested symmetry plane (x or y) by rotating STL (rollPitchYaw in degrees).
    case "$sym_plane_axis" in
      x) cyan "Rotating semi-cylinder STL: flat y=0 -> x=0"; surfaceTransformPoints -rollPitchYaw "(0 -90 0)" boundary.stl boundary.rot.stl; mv -f boundary.rot.stl boundary.stl ;;
      y) : ;; # already flat at y=0
    esac

    cyan "Combining STL..."
    cat "$STL" boundary.stl > "$COMB"
    ;;

  *)
    red "Unknown enclosure '$domain'"; exit 1;;
esac

# 4) Feature edges on combined STL
angle="$(ask "surfaceFeatureEdges angle (deg)" "90")"; ensure_number_strict "$angle"
FMS="${BASE}.fms"
cyan "Extracting feature edges from $COMB ..."
surfaceFeatureEdges -angle "$angle" "$COMB" "$FMS"
green "Feature file: $FMS"

# 5) Build patch list STRICTLY from combined STL
mapfile -t patches_unique < <(
  grep -i '^solid ' "$COMB" \
  | awk '{sub(/\r$/,"");$1="";gsub(/^ +/,"");print}' \
  | sed 's/ /_/g; s/[^A-Za-z0-9_]/_/g; s/^_//; s/_$//' \
  | awk 'length>0' \
  | sort -u || true
)

# 6) meshDict
mkdir -p system; MESHDICT="system/meshDict"; tmp="$(mktemp)"; exec 3>"$tmp"

cat >&3 <<'HDR'
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                |
| \\      /  F ield         | cfMesh: A library for mesh generation          |
|  \\    /   O peration     |                                                |
|   \\  /    A nd           |                                               |
|    \\/     M anipulation  |                                               |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version   2.0;
    format    ascii;
    class     dictionary;
    location  "system";
    object    meshDict;
}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

HDR
printf 'surfaceFile "%s";\n\n' "$FMS" >&3

# 7) max/min cell
maxCell="$(ask "maxCellSize (or 'not')" "0.25")"; ensure_number "$maxCell"
minCell="$(ask "minCellSize (or 'not')" "0.25")"; ensure_number "$minCell"
[[ "$maxCell" != "not" ]] && printf 'maxCellSize %s; //[m]\n' "$maxCell" >&3
[[ "$minCell" != "not" ]] && printf 'minCellSize %s; //[m]\n' "$minCell" >&3

# 8) boundaryCellSize
ans="$(ask "Use boundaryCellSize? (y/n)" "n")"
case "$ans" in
  [Yy]*)
    bth="$(ask "boundaryCellSizeRefinementThickness (or 'not')" "0.05")"; ensure_number "$bth"
    bcs="$(ask "boundaryCellSize (or 'not')" "0.25")"; ensure_number "$bcs"
    [[ "$bth" != "not" ]] && printf 'boundaryCellSizeRefinementThickness %s; //[m]\n' "$bth" >&3
    [[ "$bcs" != "not" ]] && printf 'boundaryCellSize %s; //[m]\n' "$bcs" >&3
    echo >>"$tmp"
  ;;
esac

# 9) localRefinement
ans="$(ask "localRefinement? (y/n)" "n")"
case "$ans" in
  [Yy]*)
    echo "localRefinement" >&3; echo "{" >&3; echo >>"$tmp"
    while true; do
      list_patches
      sel="$(ask "Enter patch name or 'getout' (or 'parapatch')" "getout")"
      maybe_run_parapatch "$sel" && continue
      [[ "$sel" == "getout" ]] && break
      yn="$(ask "cellSize? (y/n)" "n")"
      case "$yn" in
        [Yy]*)
          cs="$(ask "cellSize value (or 'not')" "0.01")"; ensure_number "$cs"
          [[ "$cs" == "not" ]] && continue
          cat >&3 <<EOF
    $sel
    {
        cellSize $cs;
    }

EOF
        ;;
        *)
          arl="$(ask "additionalRefinementLevels (or 'not')" "3")"; ensure_number "$arl"
          rth="$(ask "refinementThickness [m] (or 'not')" "0.05")"; ensure_number "$rth"
          [[ "$arl" == "not" && "$rth" == "not" ]] && continue
          printf "    %s\n    {\n" "$sel" >&3
          [[ "$arl" != "not" ]] && printf "        additionalRefinementLevels %s;\n" "$arl" >&3
          [[ "$rth" != "not" ]] && printf "        refinementThickness %s; //[m]\n" "$rth" >&3
          printf "    }\n\n" >&3
        ;;
      esac
    done
    echo "}" >&3; echo >>"$tmp"
  ;;
esac

# 10) objectRefinements
ans="$(ask "objectRefinements? (y/n)" "n")"
case "$ans" in
  [Yy]*)
    echo "objectRefinements" >&3; echo "{" >&3
    while true; do
      echo "Shapes: box, sphere, line, cone. Type 'getout' to finish."
      shape="$(ask "Shape" "getout")"; [[ "$shape" == "getout" ]] && break
      name="$(ask "Name" "${shape}1")"
      case "$shape" in
        box)
          cx="$(ask "centre x (or 'not')" "0")"; cy="$(ask "centre y (or 'not')" "0")"; cz="$(ask "centre z (or 'not')" "0")"
          lx="$(ask "lengthX (or 'not')" "1")"; ly="$(ask "lengthY (or 'not')" "1")"; lz="$(ask "lengthZ (or 'not')" "1")"
          cs="$(ask "cellSize (or 'not')" "0.01")"
          ensure_number "$cx"; ensure_number "$cy"; ensure_number "$cz"; ensure_number "$lx"; ensure_number "$ly"; ensure_number "$lz"; ensure_number "$cs"
          if [[ "$cx" == "not" && "$cy" == "not" && "$cz" == "not" && "$lx" == "not" && "$ly" == "not" && "$lz" == "not" && "$cs" == "not" ]]; then continue; fi
          echo ""; echo "    $name" >&3; echo "    {" >&3; echo "        type    box;" >&3
          if [[ "$cx" != "not" || "$cy" != "not" || "$cz" != "not" ]]; then echo "        centre  ( ${cx/not/0} ${cy/not/0} ${cz/not/0} );" >&3; fi
          [[ "$lx" != "not" ]] && echo "        lengthX $lx;" >&3
          [[ "$ly" != "not" ]] && echo "        lengthY $ly;" >&3
          [[ "$lz" != "not" ]] && echo "        lengthZ $lz;" >&3
          [[ "$cs" != "not" ]] && echo "        cellSize $cs;" >&3
          echo "    }" >&3
        ;;
        sphere)
          cx="$(ask "centre x (or 'not')" "0")"; cy="$(ask "centre y (or 'not')" "0")"; cz="$(ask "centre z (or 'not')" "0")"
          r="$(ask "radius (or 'not')" "0.1")"; cs="$(ask "cellSize (or 'not')" "0.01")"; rt="$(ask "refinementThickness (or 'not')" "not")"
          ensure_number "$cx"; ensure_number "$cy"; ensure_number "$cz"; ensure_number "$r"; ensure_number "$cs"; ensure_number "$rt"
          if [[ "$cx" == "not" && "$cy" == "not" && "$cz" == "not" && "$r" == "not" && "$cs" == "not" && "$rt" == "not" ]]; then continue; fi
          echo ""; echo "    $name" >&3; echo "    {" >&3; echo "        type   sphere;" >&3
          if [[ "$cx" != "not" || "$cy" != "not" || "$cz" != "not" ]]; then echo "        centre ( ${cx/not/0} ${cy/not/0} ${cz/not/0} );" >&3; fi
          [[ "$r"  != "not" ]] && echo "        radius $r;" >&3
          [[ "$cs" != "not" ]] && echo "        cellSize $cs;" >&3
          [[ "$rt" != "not" ]] && echo "        refinementThickness $rt;" >&3
          echo "    }" >&3
        ;;
        line)
          x0="$(ask "p0 x (or 'not')" "-0.5")"; y0="$(ask "p0 y (or 'not')" "0")"; z0="$(ask "p0 z (or 'not')" "0")"
          x1="$(ask "p1 x (or 'not')" "0.5")"; y1="$(ask "p1 y (or 'not')" "0")"; z1="$(ask "p1 z (or 'not')" "0")"
          rad="$(ask "radius (or 'not')" "0.05")"; cs="$(ask "cellSize (or 'not')" "0.01")"; rt="$(ask "refinementThickness (or 'not')" "not")"
          ensure_number "$x0"; ensure_number "$y0"; ensure_number "$z0"; ensure_number "$x1"; ensure_number "$y1"; ensure_number "$z1"; ensure_number "$rad"; ensure_number "$cs"; ensure_number "$rt"
          if [[ "$x0" == "not" && "$y0" == "not" && "$z0" == "not" && "$x1" == "not" && "$y1" == "not" && "$z1" == "not" && "$rad" == "not" && "$cs" == "not" && "$rt"  == "not" ]]; then continue; fi
          echo ""; echo "    $name" >&3; echo "    {" >&3; echo "        type   line;" >&3
          if [[ "$x0" != "not" || "$y0" != "not" || "$z0" != "not" ]]; then echo "        p0     ( ${x0/not/0} ${y0/not/0} ${z0/not/0} );" >&3; fi
          if [[ "$x1" != "not" || "$y1" != "not" || "$z1" != "not" ]]; then echo "        p1     ( ${x1/not/0} ${y1/not/0} ${z1/not/0} );" >&3; fi
          [[ "$rad" != "not" ]] && echo "        radius $rad;" >&3
          [[ "$cs"  != "not" ]] && echo "        cellSize $cs;" >&3
          [[ "$rt"  != "not" ]] && echo "        refinementThickness $rt;" >&3
          echo "    }" >&3
        ;;
        cone)
          x0="$(ask "p0 x (or 'not')" "0")"; y0="$(ask "p0 y (or 'not')" "0")"; z0="$(ask "p0 z (or 'not')" "0")"
          x1="$(ask "p1 x (or 'not')" "0")"; y1="$(ask "p1 y (or 'not')" "0")"; z1="$(ask "p1 z (or 'not')" "0.8")"
          r0="$(ask "radius0 (or 'not')" "0.15")"; r1="$(ask "radius1 (or 'not')" "0.03")"; cs="$(ask "cellSize (or 'not')" "0.01")"
          ensure_number "$x0"; ensure_number "$y0"; ensure_number "$z0"; ensure_number "$x1"; ensure_number "$y1"; ensure_number "$z1"; ensure_number "$r0"; ensure_number "$r1"; ensure_number "$cs"
          if [[ "$x0" == "not" && "$y0" == "not" && "$z0" == "not" && "$x1" == "not" && "$y1" == "not" && "$z1" == "not" && "$r0" == "not" && "$r1" == "not" && "$cs" == "not" ]]; then continue; fi
          echo ""; echo "    $name" >&3; echo "    {" >&3; echo "        type    cone;" >&3
          if [[ "$x0" != "not" || "$y0" != "not" || "$z0" != "not" ]]; then echo "        p0      ( ${x0/not/0} ${y0/not/0} ${z0/not/0} );" >&3; fi
          if [[ "$x1" != "not" || "$y1" != "not" || "$z1" != "not" ]]; then echo "        p1      ( ${x1/not/0} ${y1/not/0} ${z1/not/0} );" >&3; fi
          [[ "$r0" != "not" ]] && echo "        radius0 $r0;" >&3
          [[ "$r1" != "not" ]] && echo "        radius1 $r1;" >&3
          [[ "$cs" != "not" ]] && echo "        cellSize $cs;" >&3
          echo "    }" >&3
        ;;
        *) red "Invalid shape." ;;
      esac
    done
    echo "}" >&3; echo >>"$tmp"
  ;;
esac

# 11) boundaryLayers
ans="$(ask "boundaryLayers? (y/n)" "n")"
case "$ans" in
  [Yy]*)
    echo "boundaryLayers" >&3; echo "{" >&3; echo >>"$tmp"
    echo "    patchBoundaryLayers" >&3; echo "    {" >&3; echo >>"$tmp"
    while true; do
      list_patches
      sel="$(ask "Enter patch name or 'getout' (or 'parapatch')" "getout")"
      maybe_run_parapatch "$sel" && continue
      [[ "$sel" == "getout" ]] && break
      nl="$(ask "nLayers (or 'not')" "5")"; ensure_number "$nl"
      tr="$(ask "thicknessRatio (or 'not')" "1.2")"; ensure_number "$tr"
      mflt="$(ask "maxFirstLayerThickness (or 'not')" "not")"; ensure_number "$mflt"
      ad="$(ask "allowDiscontinuity (1/0 or 'not')" "1")"; ensure_number "$ad"
      [[ "$nl" == "not" && "$tr" == "not" && "$mflt" == "not" && "$ad" == "not" ]] && continue
      printf "        %s\n        {\n" "$sel" >&3
      [[ "$nl"  != "not" ]] && printf "            nLayers %s;\n" "$nl" >&3
      [[ "$tr"  != "not" ]] && printf "            thicknessRatio %s;\n" "$tr" >&3
      [[ "$mflt" != "not" ]] && printf "            maxFirstLayerThickness %s;\n" "$mflt" >&3
      [[ "$ad"  != "not" ]] && printf "            allowDiscontinuity %s;\n" "$ad" >&3
      printf "        }\n\n" >&3
    done
    echo "    }" >&3

    ans2="$(ask "optimiseLayer? (y/n)" "y")"
    case "$ans2" in
      [Yy]*)
        echo >>"$tmp"; echo "    optimiseLayer 1;" >&3; echo >>"$tmp"
        ans3="$(ask "optimisationParameters? (y/n)" "n")"
        case "$ans3" in
          [Yy]*)
            nsn="$(ask "nSmoothNormals (or 'not')" "5")"; ensure_number "$nsn"
            rtt="$(ask "relThicknessTol (or 'not')" "0.4")"; ensure_number "$rtt"
            fsf="$(ask "featureSizeFactor (or 'not')" "0.8")"; ensure_number "$fsf"
            rcn="$(ask "reCalculateNormals (or 'not')" "1")"; ensure_number "$rcn"
            mni="$(ask "maxNumIterations (or 'not')" "5")"; ensure_number "$mni"
            if [[ "$nsn" != "not" || "$rtt" != "not" || "$fsf" != "not" || "$rcn" != "not" || "$mni" != "not" ]]; then
              echo "    optimisationParameters" >&3
              echo "    {" >&3
              [[ "$nsn" != "not" ]] && printf "        nSmoothNormals \t\t%s;\n" "$nsn" >&3
              [[ "$rtt" != "not" ]] && printf "        relThicknessTol \t\t%s;\n" "$rtt" >&3
              [[ "$fsf" != "not" ]] && printf "        featureSizeFactor \t\t%s;\n" "$fsf" >&3
              [[ "$rcn" != "not" ]] && printf "        reCalculateNormals \t\t%s;\n" "$rcn" >&3
              [[ "$mni" != "not" ]] && printf "        maxNumIterations \t\t%s;\n" "$mni" >&3
              echo "    }" >&3
            fi
          ;;
        esac
      ;;
    esac

    echo "}" >&3; echo >>"$tmp"
  ;;
esac

# 12) renameBoundary (auto symmetry mapping if symmetry exists)
ans="$(ask "renameBoundary? (y/n)" "n")"
case "$ans" in
  [Yy]*)
    echo "renameBoundary" >&3; echo "{" >&3
    defName="$(ask "defaultName (or 'not')" "defaultPatch")"
    defType="$(ask "defaultType (or 'not')" "wall")"
    [[ "$defName" != "not" ]] && printf "    defaultName     %s;\n" "$defName" >&3
    [[ "$defType" != "not" ]] && printf "    defaultType     %s;\n\n" "$defType" >&3

    echo "    newPatchNames" >&3; echo "    {" >&3
    if [[ -n "$sym_plane_axis" ]]; then
      while true; do
        list_patches
        sp="$(ask "Which patch is the symmetry plane? (or 'parapatch' / 'skip')" "symmetry")"
        maybe_run_parapatch "$sp" && continue
        [[ "$sp" == "skip" ]] && break
        if [[ -n "$sp" && "$sp" != "not" ]]; then
          cat >&3 <<EOF
        "$sp"
        {
            newName symmetry;
            type    symmetryPlane;
        }
EOF
        fi
        break
      done
    fi

    while true; do
      list_patches
      old="$(ask "Existing patch name or 'getout' (or 'parapatch')" "getout")"
      maybe_run_parapatch "$old" && continue
      [[ "$old" == "getout" ]] && break
      new="$(ask "newName (or 'not' to skip this mapping)" "${old}_new")"
      typ="$(ask "type (patch/wall/symmetryPlane/inlet/outlet etc., or 'not')" "patch")"
      [[ "$new" == "not" && "$typ" == "not" ]] && continue
      echo "        \"$old\"" >&3; echo "        {" >&3
      [[ "$new" != "not" ]] && printf "            newName %s;\n" "$new" >&3
      [[ "$typ" != "not" ]] && printf "            type    %s;\n" "$typ" >&3
      echo "        }" >&3
    done
    echo "    }" >&3; echo "}" >&3; echo >>"$tmp"
  ;;
esac

# 13) Footer
cat >&3 <<'FTR'
// ************************************************************************* //
FTR

exec 3>&-
mv "$tmp" "$MESHDICT"
green "meshDict written: $MESHDICT"
bold "Done."
