// FEM çekirdeği ANALİTİK kıyaslarla doğrulanır (CAE test disiplini: analitik
// çözümü olmayan sonuca güvenilmez).
#include <cmath>
#include <cstdio>
#include <vector>

#include "fiko/cae/fem_hex.hpp"
#include "fiko_test.hpp"

using namespace fiko;
using namespace fiko::cae;

namespace {

HexKafes cubuk(u32 nx, u32 ny, u32 nz, f64 h) {
  HexKafes kf;
  kf.nx = nx;
  kf.ny = ny;
  kf.nz = nz;
  kf.h = h;
  kf.dolu.assign(nx * ny * nz, 1);
  return kf;
}

}  // namespace

// Eleman rijitliği: simetrik olmalı ve RİJİT CİSİM ötelemesinde sıfır kuvvet
// üretmeli (satır toplamları = 0). Bu iki koşul bozuksa tüm çözüm çöp olur.
FIKO_TEST(fem_eleman_rijitlik_simetrik_ve_rijit_cisim_sifir) {
  const std::vector<f64> ke = hexElemanRijitlik(2.0, {210000.0, 0.3});
  FIKO_KONTROL(ke.size() == 576u);
  f64 enBuyukAsimetri = 0.0;
  for (int a = 0; a < 24; ++a)
    for (int b = 0; b < 24; ++b)
      enBuyukAsimetri = std::max(enBuyukAsimetri, std::abs(ke[a * 24 + b] - ke[b * 24 + a]));
  FIKO_KONTROL(enBuyukAsimetri < 1e-6);

  f64 enBuyukSatir = 0.0;
  for (int a = 0; a < 24; ++a) {
    f64 t = 0.0;
    for (int n = 0; n < 8; ++n) t += ke[a * 24 + 3 * n + (a % 3)];  // aynı yönde öteleme
    enBuyukSatir = std::max(enBuyukSatir, std::abs(t));
  }
  FIKO_KONTROL(enBuyukSatir < 1e-6);

  f64 kosegenEnKucuk = 1e300;
  for (int a = 0; a < 24; ++a) kosegenEnKucuk = std::min(kosegenEnKucuk, ke[a * 24 + a]);
  FIKO_KONTROL(kosegenEnKucuk > 0.0);  // pozitif tanımlılık işareti
}

// ANALİTİK 1: eksenel çekme. σ = F/A TAM çıkmalı (sabit gerilme alanı hex8 ile
// birebir temsil edilir). Yer değiştirme δ = F·L/(A·E) analitiğinden BİRAZ KÜÇÜK
// çıkar: kök yüzü 3 YÖNDE ankastre tutulduğu için Poisson daralması engellenir ve
// kök bölgesi yerel olarak sertleşir (ölçüldü 07-21: L=20h'de %1.12).
// Bu YEREL bir etkidir → boy uzadıkça bağıl hata ORANTILI küçülmeli. Test bunu
// ölçer: sapma sabit bir "yanlış rijitlik" olsaydı boyla değişmezdi.
namespace {

struct CekmeOlcum {
  f64 sapma = 0.0;      // (hesap − analitik)/analitik
  f64 gerilmeHata = 0.0;
  bool yakinsadi = false;
};

CekmeOlcum eksenelCekme(u32 nx) {
  const f64 h = 2.0;
  const u32 ny = 2, nz = 2;
  const HexKafes kf = cubuk(nx, ny, nz, h);
  const ElastikMalzeme mlz{210000.0, 0.3};
  const f64 L = nx * h, A = (ny * h) * (nz * h), F = 1000.0;

  std::vector<u32> sabit;
  for (u32 k = 0; k <= nz; ++k)
    for (u32 j = 0; j <= ny; ++j) sabit.push_back(kf.dugumIndeks(0, j, k));

  // Uç yüzde ALAN PAYLI dağıtım (köşe 1/4, kenar 1/2, iç 1) → düzgün çekme.
  std::vector<DugumKuvvet> yuk;
  f64 toplam = 0.0;
  for (u32 k = 0; k <= nz; ++k)
    for (u32 j = 0; j <= ny; ++j)
      toplam += ((j == 0 || j == ny) ? 0.5 : 1.0) * ((k == 0 || k == nz) ? 0.5 : 1.0);
  for (u32 k = 0; k <= nz; ++k)
    for (u32 j = 0; j <= ny; ++j) {
      const f64 a = ((j == 0 || j == ny) ? 0.5 : 1.0) * ((k == 0 || k == nz) ? 0.5 : 1.0);
      yuk.push_back({kf.dugumIndeks(nx, j, k), {F * a / toplam, 0.0, 0.0}});
    }

  const FemSonuc s = femCoz(kf, mlz, sabit, yuk);
  CekmeOlcum o;
  o.yakinsadi = s.yakinsadi;
  if (!s.yakinsadi) return o;
  const f64 analitik = F * L / (A * mlz.elastisite);
  o.sapma = (s.yer[kf.dugumIndeks(nx, 1, 1)].x - analitik) / analitik;
  o.gerilmeHata =
      std::abs(s.vonMises[kf.hucreIndeks(nx / 2, ny / 2, nz / 2)] - F / A) / (F / A);
  std::printf("    OLCUM nx=%2u: uzama sapma=%%%.3f  gerilme hata=%%%.4f\n", nx,
              100.0 * o.sapma, 100.0 * o.gerilmeHata);
  return o;
}

}  // namespace

FIKO_TEST(fem_eksenel_cekme_analitik) {
  const CekmeOlcum kisa = eksenelCekme(10);
  const CekmeOlcum uzun = eksenelCekme(40);
  FIKO_KONTROL(kisa.yakinsadi && uzun.yakinsadi);

  // Gerilme sabit alan → TAM (yalnız çözücü toleransı kadar hata).
  FIKO_KONTROL(kisa.gerilmeHata < 1e-4);
  FIKO_KONTROL(uzun.gerilmeHata < 1e-4);

  // Yer değiştirme analitikten AZ (daha rijit) ve mutlak sapma küçük.
  FIKO_KONTROL(kisa.sapma < 0.0 && std::abs(kisa.sapma) < 0.02);
  // 4 kat boy → kök etkisi 4 kata yakın seyrelir (yerel etki kanıtı).
  const f64 oran = std::abs(kisa.sapma) / std::abs(uzun.sapma);
  std::printf("    OLCUM sapma orani (nx10/nx40) = %.2f  [~4 beklenir]\n", oran);
  FIKO_KONTROL(oran > 3.0 && oran < 5.0);
}

// ANALİTİK 2: yükün İKİ KATI → yer değiştirme de tam iki kat (lineerlik).
FIKO_TEST(fem_lineerlik_yuk_iki_kat) {
  const HexKafes kf = cubuk(6, 2, 2, 2.0);
  const ElastikMalzeme mlz{210000.0, 0.3};
  std::vector<u32> sabit;
  for (u32 k = 0; k <= 2; ++k)
    for (u32 j = 0; j <= 2; ++j) sabit.push_back(kf.dugumIndeks(0, j, k));
  const u32 uc = kf.dugumIndeks(6, 1, 1);

  const FemSonuc a = femCoz(kf, mlz, sabit, {{uc, {0.0, 0.0, -500.0}}});
  const FemSonuc b = femCoz(kf, mlz, sabit, {{uc, {0.0, 0.0, -1000.0}}});
  FIKO_KONTROL(a.yakinsadi && b.yakinsadi);
  FIKO_KONTROL(std::abs(a.yer[uc].z) > 1e-9);
  FIKO_YAKIN(b.yer[uc].z / a.yer[uc].z, 2.0, 1e-4);
}

// Kısıtsız model tekildir: çözücü "yakınsadı" DEMEMELİ (sessiz çöp sonuç yasak).
FIKO_TEST(fem_kisit_yoksa_yakinsamaz) {
  const HexKafes kf = cubuk(2, 2, 2, 1.0);
  const std::vector<u32> kisitYok;
  const FemSonuc s = femCoz(kf, {}, kisitYok, {{0, {1.0, 0.0, 0.0}}});
  FIKO_KONTROL(!s.yakinsadi);
  FIKO_KONTROL(s.yer.empty());
}

// Boşluklu kafes: dolu=0 hücreler eleman üretmez, o bölgedeki düğümler kısıtlı
// sayılır (tekillik yok) ve çözüm yine yakınsar.
FIKO_TEST(fem_bosluklu_kafes_yakinsar) {
  HexKafes kf = cubuk(8, 2, 2, 2.0);
  for (u32 k = 0; k < 2; ++k) kf.dolu[kf.hucreIndeks(4, 1, k)] = 0;  // içeride delik
  const ElastikMalzeme mlz{210000.0, 0.3};
  std::vector<u32> sabit;
  for (u32 k = 0; k <= 2; ++k)
    for (u32 j = 0; j <= 2; ++j) sabit.push_back(kf.dugumIndeks(0, j, k));
  const FemSonuc s = femCoz(kf, mlz, sabit, {{kf.dugumIndeks(8, 1, 1), {800.0, 0.0, 0.0}}});
  FIKO_KONTROL(s.yakinsadi);
  FIKO_KONTROL(s.enBuyukYer > 0.0);
  FIKO_KONTROL(s.vonMises[kf.hucreIndeks(4, 1, 0)] == 0.0);  // boş hücrede gerilme YOK
}

// İPTAL: uzun çözümde belge kapanışı beklemesin diye CG erken çıkabilmeli.
// Erken çıkışta sonuç YAKINSAMADI olarak işaretlenmeli (yarım sonuç "doğru" sanılmasın).
FIKO_TEST(fem_iptal_erken_cikar) {
  const HexKafes kf = cubuk(20, 6, 6, 1.0);
  std::vector<u32> sabit;
  for (u32 k = 0; k <= 6; ++k)
    for (u32 j = 0; j <= 6; ++j) sabit.push_back(kf.dugumIndeks(0, j, k));
  const std::vector<DugumKuvvet> yuk{{kf.dugumIndeks(20, 3, 3), {0.0, 0.0, -100.0}}};

  const FemSonuc tam = femCoz(kf, {210000.0, 0.3}, sabit, yuk);
  FIKO_KONTROL(tam.yakinsadi);
  FIKO_KONTROL(tam.adim > 1);

  int cagri = 0;
  const FemSonuc kesik =
      femCoz(kf, {210000.0, 0.3}, sabit, yuk, [&cagri] { ++cagri; return true; });
  FIKO_KONTROL(cagri > 0);            // iptal SORULDU
  FIKO_KONTROL(!kesik.yakinsadi);     // yarım sonuç "yakınsadı" DEMEZ
  FIKO_KONTROL(kesik.adim < tam.adim);  // gerçekten erken çıktı
}