% Code intended to easen the learning curve of the FEM and TopOpt methods

clear all; close all; clc;

addpath(genpath('./Programa_00'))

disp(' ')
disp('         ***************')
disp('         **   START   **')
disp('         ***************')
disp(' ')


%% --------------------------------------------------------------------- %%
%                              ** Input **                                %
%-------------------------------------------------------------------------%

% Define mesh dimensions
Lx = 20;
Ly = 60;
% Number of elements in the mesh (in x and y)
nelx = 20;
nely = 60;

% Material properties
E = 1.0;   % Young's modulus
nu = 0.3;  % Poisson's ratio
rho = 1.0; % Material density (not used in this example)

%% --------------------------------------------------------------------- %%
%                         ** Problem set up **                            %
%-------------------------------------------------------------------------%

% Initiate mesh object
mesh = Mesh();

% Generate regulard mesh grid with quadrilateral elements
mesh = GenerateQuad4RectangleMesh(mesh,[Lx Ly],[nelx nely]);

% Add Dirichlet boundary conditions
dirichlet.point_1 = [-0.001, -0.001];
dirichlet.point_2 = [mesh.Lx*1.001, 0.001];
dirichlet.dofs = [1, 2];
dirichlet.value = 0.0;
mesh = AddDirichletBC(mesh,dirichlet);

% Add Neumann boundary conditions
neumann.point_1 = [mesh.Lx*0.4999, mesh.Ly*0.9999];
neumann.point_2 = [mesh.Lx*0.5001, mesh.Ly*1.0001];
neumann.dofs = [1];
neumann.value = 1.0;
mesh = AddNeumannBC(mesh,neumann);

% Prepare elasticity study object
fea = Elasticity(mesh);

% Add material properties
fea = AddSolidMaterial(fea, E, nu, rho);

%% --------------------------------------------------------------------- %%
%                     ** Finite Element Analysis **                       %
%-------------------------------------------------------------------------%

% Assembly of load vector
fea = AssemblePointLoads(fea);

% Assembly of stiffness matrix
densities = ones(length(fea.mesh.incidence),1);
stress_case = 1; % 1 = plane stress, 2 = plane strain.
fea = AssembleStructuralK(fea,densities,stress_case);

% Solve equilibrium equation
fea = SolveStaticStructuralFEA(fea);

% Plot structural displacements
PlotStructuralDisplacements(fea,1e-1)

disp(' ')
disp('         *************')
disp('         **   END   **')
disp('         *************')
disp(' ')
