#include "fiko/cae/fem_hex.hpp"

#include <algorithm>
#include <barrier>
#include <cmath>
#include <thread>

namespace fiko::cae {

int& femEnFazlaAdim() {
  static int a = 5000;
  return a;
}

bool& konservatifKafes() {
  static bool a = true;  // ince cidar kaybı yapısal analizde sessiz felakettir
  return a;
}

bool& femParalel() {
  // ÖLÇÜLDÜ 07-21 → VARSAYILAN KAPALI. 8-renkli paralel matvec KAZANDIRMADI:
  // 2560 hücrede 129 → 206 ms (YAVAŞLADI), 8640 hücrede 651 → 565 ms (yalnız 1.15×).
  // Sebep: bariyer + iş parçacığı kurulumu HER matvec'te ödeniyor (285 iterasyon ×
  // 8 renk). Gerçek çözüm KALICI iş parçacığı havuzu (çözümün tamamı tek paralel
  // bölgede) — o yapılana dek seri yol daha hızlı. Kod ölçüm için duruyor.
  static bool a = false;
  return a;
}

bool& femBlokOnkosul() {
  // ÖLÇÜLDÜ 07-21 → VARSAYILAN KAPALI. 3×3 blok Jacobi iterasyonu yalnız **%2**
  // düşürdü (94→90, 190→186, 285→280) ve süre kazancı gürültü içinde kaldı;
  // karşılığında düğüm başına 9 f64 fazladan bellek istiyor. Çözüm birebir aynı
  // (bağıl fark ~1e-12) — yani ön-koşullayıcı doğru, ama bu ızgarada x/y/z bağı
  // beklendiği kadar güçlü DEĞİL. Asıl kazanç çok-ızgara (multigrid) tarafında.
  static bool a = false;
  return a;
}

u32& femParalelEsik() {
  static u32 a = 2000;  // hücre; altında bariyer maliyeti kazancı yer
  return a;
}

u32& asgariCidarHucre() {
  static u32 a = 4;  // bir kesitte en az 4 eleman (eğilme için pratik alt sınır)
  return a;
}

f64& femTolerans() {
  static f64 a = 1e-8;
  return a;
}

namespace {

// Hex8 düğüm sırası (yerel ξηζ ∈ {−1,+1}) — hücre (i,j,k) için
// (i+dx, j+dy, k+dz) ile birebir eşleşir.
constexpr int kYerelDx[8] = {0, 1, 1, 0, 0, 1, 1, 0};
constexpr int kYerelDy[8] = {0, 0, 1, 1, 0, 0, 1, 1};
constexpr int kYerelDz[8] = {0, 0, 0, 0, 1, 1, 1, 1};

// 2×2×2 Gauss noktaları (tam integrasyon — düşük integrasyon hourglass üretir).
constexpr f64 kG = 0.5773502691896257;  // 1/√3

// Şekil fonksiyonu türevleri (doğal koordinatta), düğüm n için.
inline void sekilTurev(f64 ksi, f64 eta, f64 zeta, int n, f64& dNdKsi, f64& dNdEta,
                       f64& dNdZeta) {
  const f64 s = kYerelDx[n] ? 1.0 : -1.0;
  const f64 t = kYerelDy[n] ? 1.0 : -1.0;
  const f64 u = kYerelDz[n] ? 1.0 : -1.0;
  dNdKsi = 0.125 * s * (1.0 + t * eta) * (1.0 + u * zeta);
  dNdEta = 0.125 * (1.0 + s * ksi) * t * (1.0 + u * zeta);
  dNdZeta = 0.125 * (1.0 + s * ksi) * (1.0 + t * eta) * u;
}

// İzotropik 6×6 malzeme matrisi (Voigt: xx yy zz xy yz zx).
void malzemeMatrisi(const ElastikMalzeme& m, f64 D[6][6]) {
  for (int a = 0; a < 6; ++a)
    for (int b = 0; b < 6; ++b) D[a][b] = 0.0;
  const f64 e = m.elastisite, nu = m.poisson;
  const f64 k = e / ((1.0 + nu) * (1.0 - 2.0 * nu));
  const f64 kayma = e / (2.0 * (1.0 + nu));
  D[0][0] = D[1][1] = D[2][2] = k * (1.0 - nu);
  D[0][1] = D[0][2] = D[1][0] = D[1][2] = D[2][0] = D[2][1] = k * nu;
  D[3][3] = D[4][4] = D[5][5] = kayma;
}

// Bir Gauss noktasında B matrisi (6×24). Düzgün küpte Jacobian = (h/2)·I.
void bMatrisi(f64 ksi, f64 eta, f64 zeta, f64 h, f64 B[6][24]) {
  for (int a = 0; a < 6; ++a)
    for (int b = 0; b < 24; ++b) B[a][b] = 0.0;
  const f64 tersJ = 2.0 / h;  // dN/dx = (2/h)·dN/dξ
  for (int n = 0; n < 8; ++n) {
    f64 dk, de, dz;
    sekilTurev(ksi, eta, zeta, n, dk, de, dz);
    const f64 dx = dk * tersJ, dy = de * tersJ, dzz = dz * tersJ;
    const int c = 3 * n;
    B[0][c + 0] = dx;
    B[1][c + 1] = dy;
    B[2][c + 2] = dzz;
    B[3][c + 0] = dy;  B[3][c + 1] = dx;   // γxy
    B[4][c + 1] = dzz; B[4][c + 2] = dy;   // γyz
    B[5][c + 0] = dzz; B[5][c + 2] = dx;   // γzx
  }
}

}  // namespace

std::vector<f64> hexElemanRijitlik(f64 h, const ElastikMalzeme& mlz) {
  std::vector<f64> ke(24 * 24, 0.0);
  if (h <= 0.0) return ke;
  f64 D[6][6];
  malzemeMatrisi(mlz, D);
  const f64 detJ = (h / 2.0) * (h / 2.0) * (h / 2.0);  // dV = detJ·dξdηdζ
  const f64 g[2] = {-kG, kG};
  f64 B[6][24], DB[6][24];
  for (int gi = 0; gi < 2; ++gi)
    for (int gj = 0; gj < 2; ++gj)
      for (int gk = 0; gk < 2; ++gk) {
        bMatrisi(g[gi], g[gj], g[gk], h, B);
        for (int a = 0; a < 6; ++a)
          for (int b = 0; b < 24; ++b) {
            f64 s = 0.0;
            for (int c = 0; c < 6; ++c) s += D[a][c] * B[c][b];
            DB[a][b] = s;
          }
        for (int a = 0; a < 24; ++a)
          for (int b = 0; b < 24; ++b) {
            f64 s = 0.0;
            for (int c = 0; c < 6; ++c) s += B[c][a] * DB[c][b];
            ke[a * 24 + b] += s * detJ;  // ağırlıklar 1·1·1
          }
      }
  return ke;
}

namespace {

// Eleman düğümlerinin GLOBAL indeksleri.
inline void elemanDugumleri(const HexKafes& kf, u32 i, u32 j, u32 k, u32 dg[8]) {
  for (int n = 0; n < 8; ++n)
    dg[n] = kf.dugumIndeks(i + kYerelDx[n], j + kYerelDy[n], k + kYerelDz[n]);
}

// Tek hücrenin katkısını y'ye ekler (ortak gövde — seri ve paralel yol AYNI kodu kullanır).
inline void hucreKatki(const HexKafes& kf, const std::vector<f64>& ke, const std::vector<f64>& x,
                       std::vector<f64>& y, u32 i, u32 j, u32 k) {
  if (!kf.dolu[kf.hucreIndeks(i, j, k)]) return;
  u32 dg[8];
  elemanDugumleri(kf, i, j, k, dg);
  f64 xe[24], ye[24];
  for (int n = 0; n < 8; ++n)
    for (int d = 0; d < 3; ++d) xe[3 * n + d] = x[3 * dg[n] + d];
  for (int a = 0; a < 24; ++a) {
    f64 s = 0.0;
    const f64* satir = &ke[a * 24];
    for (int b = 0; b < 24; ++b) s += satir[b] * xe[b];
    ye[a] = s;
  }
  for (int n = 0; n < 8; ++n)
    for (int d = 0; d < 3; ++d) y[3 * dg[n] + d] += ye[3 * n + d];
}

// y = K·x (matrix-free, eleman-eleman). Kısıtlı serbestlikler sıfırlanır.
//
// PARALEL: 8 RENKLİ şema. (i%2, j%2, k%2) aynı olan hücreler ORTAK DÜĞÜM PAYLAŞMAZ
// → aynı renk içinde eşzamanlı `y[...] +=` yarışsızdır. Renkler SIRAYLA işlenir, yani
// toplama sırası iş parçacığı SAYISINDAN bağımsız ve tekrarlanabilirdir.
// UYARI: sıra seri yoldan farklıdır → kayan nokta sonuçları son bitlerde (~1e-16 bağıl)
// değişebilir. Ayar `cae.fem_paralel` ile kapatılıp seri yola dönülebilir.
void carp(const HexKafes& kf, const std::vector<f64>& ke, const std::vector<u8>& kisitli,
          const std::vector<f64>& x, std::vector<f64>& y, u32 isSayi) {
  std::fill(y.begin(), y.end(), 0.0);
  if (isSayi <= 1) {
    for (u32 k = 0; k < kf.nz; ++k)
      for (u32 j = 0; j < kf.ny; ++j)
        for (u32 i = 0; i < kf.nx; ++i) hucreKatki(kf, ke, x, y, i, j, k);
  } else {
    std::barrier engel((std::ptrdiff_t)isSayi);
    const auto is = [&](u32 t) {
      for (int renk = 0; renk < 8; ++renk) {
        const u32 rk = (u32)(renk & 1), rj = (u32)((renk >> 1) & 1), ri = (u32)((renk >> 2) & 1);
        for (u32 k = rk + 2 * t; k < kf.nz; k += 2 * isSayi)
          for (u32 j = rj; j < kf.ny; j += 2)
            for (u32 i = ri; i < kf.nx; i += 2) hucreKatki(kf, ke, x, y, i, j, k);
        engel.arrive_and_wait();
      }
    };
    std::vector<std::thread> isler;
    isler.reserve(isSayi - 1);
    for (u32 t = 1; t < isSayi; ++t) isler.emplace_back(is, t);
    is(0);
    for (std::thread& i : isler) i.join();
  }
  for (std::size_t s = 0; s < y.size(); ++s)
    if (kisitli[s]) y[s] = 0.0;
}

}  // namespace

FemSonuc femCoz(const HexKafes& kafes, const ElastikMalzeme& mlz,
                const std::vector<DugumKisit>& kisitlar,
                const std::vector<DugumKuvvet>& yukler, const std::function<bool()>& iptal) {
  FemSonuc son;
  if (!kafes.gecerliMi() || kisitlar.empty() || yukler.empty()) return son;

  const u32 nd = kafes.dugumSayi();
  const std::size_t n = (std::size_t)nd * 3;
  const std::vector<f64> ke = hexElemanRijitlik(kafes.h, mlz);

  // Kısıt maskesi + serbestliği hiç görmeyen (bağsız) düğümler de kısıtlıdır:
  // aksi halde K tekildir ve CG sonsuza koşar.
  std::vector<u8> kisitli(n, 1);
  u32 dg[8];
  for (u32 k = 0; k < kafes.nz; ++k)
    for (u32 j = 0; j < kafes.ny; ++j)
      for (u32 i = 0; i < kafes.nx; ++i) {
        if (!kafes.dolu[kafes.hucreIndeks(i, j, k)]) continue;
        elemanDugumleri(kafes, i, j, k, dg);
        for (int m = 0; m < 8; ++m)
          for (int d = 0; d < 3; ++d) kisitli[3 * dg[m] + d] = 0;
      }
  for (const DugumKisit& kst : kisitlar) {
    if (kst.dugum >= nd) continue;
    for (int d = 0; d < 3; ++d)
      if (kst.eksenMaskesi & (1u << d)) kisitli[3 * kst.dugum + d] = 1;
  }

  std::vector<f64> b(n, 0.0);
  for (const DugumKuvvet& y : yukler) {
    if (y.dugum >= nd) continue;
    b[3 * y.dugum + 0] += y.kuvvet.x;
    b[3 * y.dugum + 1] += y.kuvvet.y;
    b[3 * y.dugum + 2] += y.kuvvet.z;
  }
  for (std::size_t s = 0; s < n; ++s)
    if (kisitli[s]) b[s] = 0.0;

  // ÖN-KOŞULLAYICI. İki seçenek:
  //  · SKALER Jacobi: M = diag(K).
  //  · BLOK Jacobi (3×3, düğüm başına): elastiklikte bir düğümün x/y/z serbestlikleri
  //    KUVVETLE bağlıdır (ν≠0 → çapraz terimler); skaler köşegen bu bağı görmez.
  //    Düğüm bloğunu tersine çevirmek iterasyonu düşürür (aşağıda ölçüldü).
  // Kısıtlı serbestlikler her iki yolda da BİRİM bırakılır.
  std::vector<f64> kosegen(n, 0.0);           // skaler yol
  std::vector<f64> blok((std::size_t)nd * 9, 0.0);  // 3×3 blok (satır-öncelikli)
  const bool blokKip = femBlokOnkosul();
  for (u32 k = 0; k < kafes.nz; ++k)
    for (u32 j = 0; j < kafes.ny; ++j)
      for (u32 i = 0; i < kafes.nx; ++i) {
        if (!kafes.dolu[kafes.hucreIndeks(i, j, k)]) continue;
        elemanDugumleri(kafes, i, j, k, dg);
        for (int m = 0; m < 8; ++m)
          for (int d = 0; d < 3; ++d) {
            const int a = 3 * m + d;
            kosegen[3 * dg[m] + d] += ke[a * 24 + a];
            if (blokKip)
              for (int e = 0; e < 3; ++e)
                blok[(std::size_t)dg[m] * 9 + d * 3 + e] += ke[a * 24 + (3 * m + e)];
          }
      }
  for (std::size_t s = 0; s < n; ++s)
    if (kisitli[s] || kosegen[s] <= 0.0) kosegen[s] = 1.0;
  if (blokKip) {
    // Blokları YERİNDE tersine çevir; kısıtlı/tekil satır → birim (skaler yola düşer).
    for (u32 d = 0; d < nd; ++d) {
      f64* b = &blok[(std::size_t)d * 9];
      const bool kisitliDugum =
          kisitli[3 * d] || kisitli[3 * d + 1] || kisitli[3 * d + 2];
      const f64 det = b[0] * (b[4] * b[8] - b[5] * b[7]) - b[1] * (b[3] * b[8] - b[5] * b[6]) +
                      b[2] * (b[3] * b[7] - b[4] * b[6]);
      if (kisitliDugum || std::abs(det) < 1e-30) {
        for (int t = 0; t < 9; ++t) b[t] = 0.0;
        b[0] = 1.0 / kosegen[3 * d];
        b[4] = 1.0 / kosegen[3 * d + 1];
        b[8] = 1.0 / kosegen[3 * d + 2];
        continue;
      }
      const f64 c[9] = {b[4] * b[8] - b[5] * b[7], b[2] * b[7] - b[1] * b[8],
                        b[1] * b[5] - b[2] * b[4], b[5] * b[6] - b[3] * b[8],
                        b[0] * b[8] - b[2] * b[6], b[2] * b[3] - b[0] * b[5],
                        b[3] * b[7] - b[4] * b[6], b[1] * b[6] - b[0] * b[7],
                        b[0] * b[4] - b[1] * b[3]};
      for (int t = 0; t < 9; ++t) b[t] = c[t] / det;
    }
  }
  // z = M⁻¹ r (tek kapı: skaler ya da blok).
  const auto onkosul = [&](const std::vector<f64>& rr, std::vector<f64>& zz) {
    if (!blokKip) {
      for (std::size_t s = 0; s < n; ++s) zz[s] = rr[s] / kosegen[s];
      return;
    }
    for (u32 d = 0; d < nd; ++d) {
      const f64* b = &blok[(std::size_t)d * 9];
      const f64 r0 = rr[3 * d], r1 = rr[3 * d + 1], r2 = rr[3 * d + 2];
      zz[3 * d] = b[0] * r0 + b[1] * r1 + b[2] * r2;
      zz[3 * d + 1] = b[3] * r0 + b[4] * r1 + b[5] * r2;
      zz[3 * d + 2] = b[6] * r0 + b[7] * r1 + b[8] * r2;
    }
  };

  // Ön-koşullu CG (u0 = 0 → r0 = b).
  std::vector<f64> u(n, 0.0), r = b, z(n), p(n), q(n);
  f64 bNorm = 0.0;
  for (f64 v : b) bNorm += v * v;
  bNorm = std::sqrt(bNorm);
  if (bNorm <= 0.0) return son;  // yük tümüyle kısıtlı düğümlerdeydi

  onkosul(r, z);
  p = z;
  f64 rz = 0.0;
  for (std::size_t s = 0; s < n; ++s) rz += r[s] * z[s];

  // Paralel eşiği: küçük modelde iş parçacığı + bariyer maliyeti kazançtan pahalı.
  const u32 isSayi = (femParalel() && kafes.hucreSayi() >= femParalelEsik())
                         ? std::max(1u, std::min(kafes.nz / 2,
                                                 std::thread::hardware_concurrency()))
                         : 1u;
  const int enFazla = femEnFazlaAdim();
  const f64 tol = femTolerans();
  int adim = 0;
  f64 bagilArtik = 1.0;
  for (; adim < enFazla; ++adim) {
    // İPTAL: 32 adımda bir sorulur (atomik okuma sıcak yolda gürültü yapmasın).
    // Erken çıkışta `yakinsadi` FALSE kalır → çağıran yarım sonucu doğru sanmaz.
    if (iptal && (adim & 31) == 0 && iptal()) break;
    carp(kafes, ke, kisitli, p, q, isSayi);
    f64 pq = 0.0;
    for (std::size_t s = 0; s < n; ++s) pq += p[s] * q[s];
    if (pq <= 0.0) break;  // K yarı-tanımlı değil (kısıt eksik) → dur
    const f64 alfa = rz / pq;
    for (std::size_t s = 0; s < n; ++s) {
      u[s] += alfa * p[s];
      r[s] -= alfa * q[s];
    }
    f64 rNorm = 0.0;
    for (f64 v : r) rNorm += v * v;
    bagilArtik = std::sqrt(rNorm) / bNorm;
    if (bagilArtik < tol) {
      ++adim;
      break;
    }
    onkosul(r, z);
    f64 rzYeni = 0.0;
    for (std::size_t s = 0; s < n; ++s) rzYeni += r[s] * z[s];
    const f64 beta = rzYeni / rz;
    rz = rzYeni;
    for (std::size_t s = 0; s < n; ++s) p[s] = z[s] + beta * p[s];
  }

  son.adim = adim;
  son.artik = bagilArtik;
  son.yakinsadi = bagilArtik < tol;
  son.yer.resize(nd);
  for (u32 s = 0; s < nd; ++s) {
    son.yer[s] = {u[3 * s + 0], u[3 * s + 1], u[3 * s + 2]};
    son.enBuyukYer = std::max(son.enBuyukYer, son.yer[s].uzunluk());
  }

  // Gerilme örnekleme noktası: ELEMAN KÖŞELERİ (ξ,η,ζ = ±1), en büyüğü alınır.
  // Merkez (0,0,0) örneklemesi eğilmede maksimumu SİSTEMATİK OLARAK KAÇIRIR: en
  // büyük lif gerilmesi yüzeydedir, merkez ise h/2 içeridedir → hata ≈ h/(2c)
  // (ölçüldü 07-21, 100×20×20 konsol: merkezle %-33 … %-15; köşeyle bu kapanır).
  f64 D[6][6];
  malzemeMatrisi(mlz, D);
  f64 Bk[8][6][24];
  for (int n = 0; n < 8; ++n)
    bMatrisi(kYerelDx[n] ? 1.0 : -1.0, kYerelDy[n] ? 1.0 : -1.0, kYerelDz[n] ? 1.0 : -1.0,
             kafes.h, Bk[n]);
  son.vonMises.assign(kafes.hucreSayi(), 0.0);
  for (u32 k = 0; k < kafes.nz; ++k)
    for (u32 j = 0; j < kafes.ny; ++j)
      for (u32 i = 0; i < kafes.nx; ++i) {
        const u32 h = kafes.hucreIndeks(i, j, k);
        if (!kafes.dolu[h]) continue;
        elemanDugumleri(kafes, i, j, k, dg);
        f64 ue[24];
        for (int m = 0; m < 8; ++m)
          for (int d = 0; d < 3; ++d) ue[3 * m + d] = u[3 * dg[m] + d];
        f64 enBuyuk = 0.0;
        for (int n = 0; n < 8; ++n) {
          f64 eps[6] = {0, 0, 0, 0, 0, 0}, sig[6] = {0, 0, 0, 0, 0, 0};
          for (int a = 0; a < 6; ++a)
            for (int b = 0; b < 24; ++b) eps[a] += Bk[n][a][b] * ue[b];
          for (int a = 0; a < 6; ++a)
            for (int c = 0; c < 6; ++c) sig[a] += D[a][c] * eps[c];
          const f64 vm = std::sqrt(
              0.5 * ((sig[0] - sig[1]) * (sig[0] - sig[1]) +
                     (sig[1] - sig[2]) * (sig[1] - sig[2]) +
                     (sig[2] - sig[0]) * (sig[2] - sig[0])) +
              3.0 * (sig[3] * sig[3] + sig[4] * sig[4] + sig[5] * sig[5]));
          enBuyuk = std::max(enBuyuk, vm);
        }
        son.vonMises[h] = enBuyuk;
        son.enBuyukVonMises = std::max(son.enBuyukVonMises, enBuyuk);
      }
  return son;
}

}  // namespace fiko::cae