% 1) PHYSICAL DATA OF THE PROBLEM
clear problem_data
% Physical domain defined as a NURBS map. We start by giving the knots,
%  the (weighted) control points and the weights for the circle of
%  radius 1.
knots = {[0 0 0 1 1 1], [0 0 0 1 1 1]};
coefs = zeros (4, 3, 3);
coefs(:,1,1) = [1, 0, 0, 1];
coefs(:,2,1) = [1/sqrt(2), -1/sqrt(2), 0, 1/sqrt(2)];
coefs(:,3,1) = [0, -1, 0, 1];
coefs(:,1,2) = [1/sqrt(2), 1/sqrt(2), 0, 1/sqrt(2)];
coefs(:,2,2) = [0, 0, 0, sqrt(2)-1];
coefs(:,3,2) = [-1/sqrt(2), -1/sqrt(2), 0, 1/sqrt(2)];
coefs(:,1,3) = [0, 1, 0, 1];
coefs(:,2,3) = [-1/sqrt(2), 1/sqrt(2), 0, 1/sqrt(2)];
coefs(:,3,3) = [-1, 0, 0, 1];
srf = nrbmak (coefs, knots);
% Then we apply a linear transform to get the circle of radius 2.
%  See the help of nrbtform to learn more.
srf = nrbtform (srf, vecscale ([2 2 0]));
% Finally, we extrude the circle to create the thin plate, and save it
%  in a .txt file that can be read by GeoPDEs.
vol = nrbextrude (srf, [0 0 0.02]);
nrbexport (vol, 'geo_thin_circular_plate.txt')
problem_data.geo_name = 'geo_thin_circular_plate.txt';

% Boundary conditions. We clamp the lateral sides (1 to 4), and leave free
%  the top and bottom sides (6 and 5, respectively).
problem_data.nmnn_sides   = [5 6];
problem_data.drchlt_sides = [1 2 3 4];
problem_data.symm_sides   = [];

% We give the same physical parameters as in Cottrell et al., 2006.
%  Notice in particular the role of the mass density.
E  =  30e6; nu = 0.2;
rho = 2.320;
problem_data.lambda_lame = @(x, y, z) ((nu*E)/((1+nu)*(1-2*nu)) * ones (size (x)));
problem_data.mu_lame = @(x, y, z) (E/(2*(1+nu)) * ones (size (x)));
problem_data.mass_density = @(x, y, z) rho * ones (size (x));

% There are no source terms, and the boundary conditions are all homogeneous.


% 2) CHOICE OF THE DISCRETIZATION PARAMETERS
%  We solve with degree 2, C^1 continuity, in a 6x6x1 mesh.
clear method_data
method_data.degree     = [2 2 2];     % Degree of the basis functions
method_data.regularity = [1 1 1];     % Regularity of the basis functions
method_data.nsub       = [6 6 1];     % Number of subdivisions
method_data.nquad      = [3 3 3];     % Points for the Gaussian quadrature rule

% 3) CALL TO THE SOLVER
%  Here we call the function that we have implemented in exercise 3 (b).
[geometry, msh, space, eigv, eigf] = solve_linear_elasticity_3d_eig (problem_data, method_data);

% 4) POST-PROCESSING.
% Exporting to Paraview is as in the other examples, but now we export each
%  eigenfunction independently, in a different file. It would be also
%  possible to write a function to export the six modes in one single file.
vtk_pts = {linspace(0, 1, 21), linspace(0, 1, 21), linspace(0, 1, 5)};
for ii = 1:numel (eigv)
  sp_to_vtk (eigf(:,ii), space, geometry, vtk_pts, ...
     sprintf ('thin_plate_mode%i_displacement.vts', ii), 'Mode')
end

disp ('Computed natural frequencies')
disp (sqrt (eigv));