%% 1D FEM Example
% We choose the 1D Helmholtz equation as a second-order differential equation to solve using Finite Element Analysis (FEM).
%
% The Helmholtz equation is the eigenvalue problem for the Laplacian.
% \nabla^2 f=-k^2 f\;
% where \nabla^2 is the laplacian 'div of grad'. In the context of waves, k is the wavenumber (periods per unit length).
% Related:
% Helmholtz' equation: \nabla^2 f=-k^2 f\;
% Particular cases of Helmholtz' equation:
% Laplace' equation: \nabla^2 f=h\left(x,y,z\right)
% Poisson's equation: \nabla^2 f=0
%
% In this example we have a transversal wave in a 1D medium (though many other phenoma may be modelled using these equations):
% -\frac{d}{\textrm{dx}}\left(\alpha \frac{\textrm{df}}{\textrm{dx}}\right)+\beta f=s
% f\left(a\right)=f_a ,f\left(b\right)=f_b
% The following parameters are prescribed:
% \alpha \left(x\right)=\frac{1}{\mu \left(x\right)}
% \beta \left(x\right)=j\omega \sigma \left(x\right)-\omega^2 \epsilon \left(x\right)
% s\left(x\right)=-j\omega J_z \left(x\right), this vanishes in the absence of current-carrying conductors.
%
%% Subdivide solution domain into elements
% In this case, we simply divide our 1D domain into equally sized cells.

a = -2
b = 5
n = 8 % Number of nodes -> n-1 elements
x = linspace(a,b,n)

%% Basis Functions
% Choose nodal basis functions \phi_i to approximate the solution at the nodes.
% f\left(r\right)\approx \sum_{i=1} f_i {\;\phi }_i \left(r\right), where f_i are unknown coefficients/weights multiplying the basis functions.
% In this case, we use "tent functions" \phi_i =1,\phi_{j\not= i} =0.
% For brevity, this can be expressed as a diagonal matrix, one row corresponding to a single basis function (there is one for each node, defined across the solution domain).

basis_fun = diag(ones(n,1));
for i=1:n
    hold on 
    plot(x, basis_fun(i,:))
end
xlabel('x[m]'); ylabel('basis func')

% In this example f\left(x\right)=\sum_{j=1}^8 f_j \phi_j \left(x\right)
% In this case, the expansions coefficients are the values of f at the nodes. ie.f\left(x_i \right)=f_i.
% Furthermore, the boundary values f_a ,f_b are known.
% Choose test, or weighting functions w_i for the unknown coefficients. Applying Galerkin's Method, we choose the weights to be the same as the basis functions w_i \left(x\right)=\phi_i \left(x\right), i=2,3...7 (endpoints excluded because corresponding function values are known).

weights = basis_fun(2:7,:) % basis functions 2 through 7.

%
% Example for importing geometry .stl file (eg. made in blender), and making a mesh.

inputSTL = ['curved_plane.stl']
model = createpde(1)
importGeometry(model, inputSTL)
generateMesh(model, 'Hmax', 0.09, 'GeometricOrder', 'linear') % 'Hmax' sets target maximum element size.
figure
pdeplot(model)

%