% the coordinate index is from 1~Nx(from left to right) for the bottom line
% and then from Nx+1~2*Nx  (from left to right) for the next line above
% and then next
xmin=0;xmax=400e3;
Nx=201;Ny=201;
x0=linspace(xmin,xmax,Nx);
ymax=200e3;
ymin=2e3*sin(x0/xmax*2*pi);k=0;
for i1=1:Nyfor i2=1:Nxk=k+1;Coord(k,:)=[x0(i2),ymin(i2)+(ymax-ymin(i2))*(i1-1)/(Ny-1)];    end
end
save coordinates.dat Coord -ascii% the element index
k=0;
elements3=zeros((Nx-1)*(Ny-1)*2,3);
for i1=1:Ny-1for i2=1:Nx-1k=k+1;ijm1=i2+(i1-1)*Nx;ijm2=i2+1+(i1-1)*Nx;ijm3=i2+1+i1*Nx;elements3(k,:)=[ijm1,ijm2,ijm3];k=k+1;ijm1=i2+1+i1*Nx;ijm2=i2+i1*Nx;ijm3=i2+(i1-1)*Nx;elements3(k,:)=[ijm1,ijm2,ijm3];end
end
%elements3=delaunay(Coord(:,1),Coord(:,2));
save elements3.dat elements3 -ascii% The direchlet boundary condition (index of the two end nodes for each boundary line)
boundary=zeros(2*(Nx-1),2);
temp1=1:Nx-1;
temp2=2:Nx;
temp3=1:Nx-1;
boundary(temp3',:)=[temp1',temp2'];
temp1=Ny*Nx:-1:Ny*Nx-Nx+2;
temp2=temp1-1;
temp3=temp3+Nx-1;
boundary(temp3',:)=[temp1',temp2'];
save dirichlet.dat boundary -ascii% The Neuemman boundary condition (index of the two end nodes for each boundary line)
boundary=zeros(2*(Ny-1),2);
temp1=Nx:Nx:Nx*(Ny-1);
temp2=temp1+Nx;
temp3=1:Nx-1;
boundary(temp3',:)=[temp1',temp2'];
temp1=1:Nx:Nx*(Ny-2)+1;
temp2=temp1+Nx;
temp3=temp3+Nx-1;
boundary(temp3',:)=[temp1',temp2'];
save neumann.dat boundary -ascii

function value = u_d ( u)%*****************************************************************************80%%% U_D evaluates the Dirichlet boundary conditions.%%%    The user must supply the appropriate routine for a given problem%%%  Parameters:%%    Input, real U(N,M), contains the M-dimensional coordinates of N points.%%    Output, VALUE(N), contains the value of the Dirichlet boundary%    condition at each point.%value = zeros ( size ( u, 1 ), 1 );value(end-length(value)/2+1:end)=1350;
temp1=length(value)/2-1;value(1:length(value)/2)=6.5e-3*2e3*sin((0:temp1)/temp1*2*pi);

function value = g ( u )%*****************************************************************************80%%% G evaluates the outward normal values assigned at Neumann boundary conditions.%%%    This routine must be changed by the user to reflect a particular problem.%%  Parameters:%%    Input, real U(N,M), contains the M-dimensional coordinates of N points.%%    Output, VALUE(N), contains the value of outward normal at each point%    where a Neumann boundary condition is applied.%value = zeros ( size ( u, 1 ), 1 );returnend

function value = f ( u )%*****************************************************************************80%%% F evaluates the right hand side of Laplace's equation. NOtice that F is qv/K instead of qv.%%%    This routine must be changed by the user to reflect a particular problem.%%%  Parameters:%%    Input, real U(N,M), contains the M-dimensional coordinates of N points.%%    Output, VALUE(N), contains the value of the right hand side of Laplace's%    equation at each of the points.%n = size ( u, 1 );value(1:n) =0;% ones(n,1);%2.0 * pi * pi * sin ( pi * u(1:n,1) ) .* sin ( pi * u(1:n,2) );

function M = stima3 ( vertices )%*****************************************************************************80
%
%% STIMA3 determines the local stiffness matrix for a triangular element.
%
%  Discussion:
%
%    Although this routine is intended for 2D usage, the same formulas
%    work for tetrahedral elements in 3D.  The spatial dimension intended
%    is determined implicitly, from the spatial dimension of the vertices.
%
%  Parameters:
%
%    Input, real VERTICES(1:(D+1),1:D), contains the D-dimensional
%    coordinates of the vertices.
%
%    Output, real M(1:(D+1),1:(D+1)), the local stiffness matrix
%    for the element.
%d = size ( vertices, 2 );D_eta = [ ones(1,d+1); vertices' ] \ [ zeros(1,d); eye(d) ];
M = det ( [ ones(1,d+1); vertices' ] ) * D_eta * D_eta' / prod ( 1:d );

%*****************************************************************************80
%
%% Aapplies the finite element method to steady heat conduction equation,
%  that is Poission's equation
%
%
%    The user supplies datafiles that specify the geometry of the region
%    and its arrangement into triangular or quadrilateral elements, and
%    the location and type of the boundary conditions, which can be any
%    mixture  of Neumann and Dirichlet.
%
%    The unknown state variable U(x,y) is assumed to satisfy
%    Poisson's equation (steady heat conduction equation):
%      -Uxx(x,y) - Uyy(x,y) = F(x,y) in Omega
%    with Dirichlet boundary conditions
%      U(x,y) = U_D(x,y) on Gamma_D
%    and Neumann boundary conditions on the outward normal derivative:
%      Un(x,y) = G(x,y) on Gamma_N
%    If Gamma designates the boundary of the region Omega,
%    then we presume that
%      Gamma = Gamma_D + Gamma_N
%    F(x,y) is qv/K,here qv means volumetric heat generation rate
%    but the user is free to determine which boundary conditions to
%    apply.
%
%    The code uses piecewise linear basis functions for triangular elements,
%    and piecewise isoparametric bilinear basis functions for quadrilateral
%    elements.
%
%
%clear%
%  Read the nodal coordinate data file.
%load coordinates.dat;
%
%  Read the triangular element data file.
%load elements3.dat;
%
%  Read the quadrilateral element data file.
%
%  load elements4.dat;
%
%  Read the Neumann boundary condition data file.
%  I THINK the purpose of the EVAL command is to create an empty NEUMANN array
%  if no Neumann file is found.
%eval ( 'load neumann.dat;', 'neumann=[];' );
%
%  Read the Dirichlet boundary condition data file.
%load dirichlet.dat;A = sparse ( size(coordinates,1), size(coordinates,1) );b = sparse ( size(coordinates,1), 1 );
%
%  Assembly.
%%{for j = 1 : size(elements3,1)A(elements3(j,:),elements3(j,:)) = A(elements3(j,:),elements3(j,:)) ...+ stima3(coordinates(elements3(j,:),:));end
%%}
%{for j = 1 : size(elements4,1)A(elements4(j,:),elements4(j,:)) = A(elements4(j,:),elements4(j,:)) ...+ stima4(coordinates(elements4(j,:),:));end
%}
%  Volume Forces.
%
% from the center of each element to Nodes
% Notice that the result of f here means qv/K instead of qv
%%{for j = 1 : size(elements3,1)b(elements3(j,:)) = b(elements3(j,:)) ...+ det( [1,1,1; coordinates(elements3(j,:),:)'] ) * ...f(sum(coordinates(elements3(j,:),:))/3)/6;end
%%}
%{for j = 1 : size(elements4,1)b(elements4(j,:)) = b(elements4(j,:)) ...+ det([1,1,1; coordinates(elements4(j,1:3),:)'] ) * ...f(sum(coordinates(elements4(j,:),:))/4)/4;end
%}
%  Neumann conditions.
%if ( ~isempty(neumann) )for j = 1 : size(neumann,1)b(neumann(j,:)) = b(neumann(j,:)) + ...norm(coordinates(neumann(j,1),:) - coordinates(neumann(j,2),:)) * ...g(sum(coordinates(neumann(j,:),:))/2)/2;endend
%
%  Determine which nodes are associated with Dirichlet conditions.
%  Assign the corresponding entries of U, and adjust the right hand side.
%u = sparse ( size(coordinates,1), 1 );BoundNodes = unique ( dirichlet );u(BoundNodes) = u_d ( coordinates(BoundNodes,:));b = b - A * u;
%
%  Compute the solution by solving A * U = B for the remaining unknown values of U.
%FreeNodes = setdiff ( 1:size(coordinates,1), BoundNodes );u(FreeNodes) = A(FreeNodes,FreeNodes) \ b(FreeNodes);
%
%  Graphic representation.
%
%  show ( elements3, elements4, coordinates, full ( u ) );figure
%%{trisurf ( elements3, coordinates(:,1)*1e-3, coordinates(:,2)*1e-3, full ( u ));
%%}
%{
trisurf ( elements4, coordinates(:,1), coordinates(:,2), u')
%}
shading interp
xlabel('x')

function value = f ( u )%*****************************************************************************80%%% F evaluates the right hand side of Laplace's equation. NOtice that F is qv/K instead of qv.%%%    This routine must be changed by the user to reflect a particular problem.%%%  Parameters:%%    Input, real U(N,M), contains the M-dimensional coordinates of N points.%%    Output, VALUE(N), contains the value of the right hand side of Laplace's%    equation at each of the points.%n = size ( u, 1 );value(1:n) =4.6080e-7;% ones(n,1);%2.0 * pi * pi * sin ( pi * u(1:n,1) ) .* sin ( pi * u(1:n,2) );