"""FEACORE Euler-Bernoulli beam benchmark. Units: N, mm, MPa.
Run with numpy and matplotlib. Generates result plots and convergence CSV.
This is a 1D linear beam FE verification example, not a 3D solid simulation.
"""
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import csv
from pathlib import Path
out=Path(__file__).parent
E=210000.;L=200.;b=20.;h=10.;I=b*h**3/12;q=-1.
x=np.linspace(0,L,401)
exact=q*x*x*(6*L*L-4*L*x+x*x)/(24*E*I)
rows=[]
fig,axs=plt.subplots(1,2,figsize=(12,4.5),layout='constrained')
for n in [1,2,4,8]:
 le=L/n;K=np.zeros((2*n+2,2*n+2));f=np.zeros(2*n+2)
 k=E*I/le**3*np.array([[12,6*le,-12,6*le],[6*le,4*le**2,-6*le,2*le**2],[-12,-6*le,12,-6*le],[6*le,2*le**2,-6*le,4*le**2]])
 for j in range(n):
  ix=np.arange(2*j,2*j+4);K[np.ix_(ix,ix)]+=k;f[ix]+=q*le/12*np.array([6,le,6,-le])
 u=np.zeros(2*n+2);u[2:]=np.linalg.solve(K[2:,2:],f[2:]);v=[]
 for xx in x:
  j=min(int(xx/le),n-1);t=(xx-j*le)/le
  N=np.array([1-3*t*t+2*t**3,le*(t-2*t*t+t**3),3*t*t-2*t**3,le*(-t*t+t**3)])
  v.append(N@u[2*j:2*j+4])
 error=max(abs(np.array(v)-exact));r=K@u-f
 rows.append([n,u[-2],error,100*error/abs(exact[-1]),r[0],r[1]])
 axs[0].plot(x,v,label=f'{n} elements',lw=1.3)
 assert abs(r[0]+q*L)<1e-6 and abs(r[1]+q*L*L/2)<1e-5
axs[0].plot(x,exact,'k--',label='Analytical',lw=1.4);axs[0].set(xlabel='Position x (mm)',ylabel='Deflection (mm)',title='Deflection: beam FE vs analytical');axs[0].legend(fontsize=8)
axs[1].loglog([r[0] for r in rows],[r[2] for r in rows],'o-',color='#164fe5');axs[1].set(xlabel='Number of beam elements',ylabel='Maximum field error (mm)',title='Convergence of the displacement field')
for a in axs:a.grid(alpha=.2)
fig.savefig(out/'beam-results.png',dpi=180)
with open(out/'beam-results.csv','w') as file:
 w=csv.writer(file);w.writerow(['elements','tip_mm','max_field_error_mm','error_pct_tip','reaction_N','reaction_moment_Nmm']);w.writerows(rows)
print(rows)