% FEM routines in MATLAB
% This .m file is to be modified accordingly
% One can use either the GUI or the following commands.
% Please check with Matlab Help.

%%%%%% Geometry generation %%%%%%

% One can generate variables gd sf ns by using pdecirc(xc,yc,radius) plus
% and exporting the data. 

[dl,bt]=decsg(gd,sf,ns);        % Constructs the required format for meshing

[p,e,t]=initmesh(dl,'hgrad',1.1,'jiggle','mean','jiggleiter',Inf);     % Initial mesh of the domain
p=jigglemesh(p,e,t,'opt','mean','iter',Inf);                           % Improves the quality of the mesh without refining by moving the nodes
[p,e,t]=refinemesh(dl,p,e,t,'regular');                                % Refines the mesh according to the criteria in ' '

numnod = length(p);             % number of nodes


figure;
pdegplot(dl);              % plots domain
figure;
pdemesh(p,e,t)             %plots mesh


%%%%%% Stiffness K and Mass M matrices %%%%%%%%

% Check on the MATLAB help for the definition of each of the coefficients:

c=zeros(1,length(t));       
a=ones(1,length(t));
f=zeros(1,length(t));

K = sparse(numnod,numnod); 
M = sparse(numnod,numnod); 

[K,M,unused]=assema(p,t,c,a,f);         % Gives the matrices


%%%%%%%% DIFFICULTY! %%%%%%%%%

% You need to implement boundary conditions. For this, you must
%   - Identify nodes on the boundary
%   - Match approximations on the line integral with the volume
%   - Construct the associated matrix B for the boundary degrees of freedom
%   - Construct an adequate right-hand side


%%%%%%% Solution %%%%%%%%%%

% Here you can use the command \ to solve the system



Z = (-K + k^2*M + B);           % Complete matrix
b                               % right-hand side
u = Z\b;                        % Solution