%*****************************************************************************
%
%    The unknown state variable U(x,y) is assumed to satisfy
%    Poisson's equation (steady heat conduction equation):
%      -K(Uxx(x,y) + Uyy(x,y)) = qv(x,y) in Omega
%    here qv means volumetric heat generation rate,K is thermocal
%    conductivity
%
%    clearclose all
%  Read the nodal coordinate data file.load coordinates.dat;
%  Read the triangular element data file.load elements3.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;
%  Read the Convective heat transfer boundary condition data file.eval ( 'load Convective.dat;', 'Convective=[];' );
alpha=100;% Convective heat transfer coefficient
tf=20;% ambient temperature
q=-2;% Surface heat flux at Neumann boundary condition
K=200;% Thermal conductivityA = sparse ( size(coordinates,1), size(coordinates,1) );b = sparse ( size(coordinates,1), 1 );
%
%  Assembly.if ( ~isempty(Convective) )for j = 1 : size(elements3,1)A(elements3(j,:),elements3(j,:)) = A(elements3(j,:),elements3(j,:)) ...+ stima3(coordinates(elements3(j,:),:),K) ...+ stima3_1(coordinates(elements3(j,:),:),elements3(j,:),Convective,alpha);endelsefor j = 1 : size(elements3,1)A(elements3(j,:),elements3(j,:)) = A(elements3(j,:),elements3(j,:)) ...+ stima3(coordinates(elements3(j,:),:),K);endend%  Volume Forces.
%
% from the center of each element to Nodes
% Notice that the result of f here means qv/K instead of qvfor 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
%  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,q)/2;endend%  Convective heat transfer boundary condition.if ( ~isempty(Convective) )for j = 1 : size(Convective,1)b(Convective(j,:)) = b(Convective(j,:)) + ...norm(coordinates(Convective(j,1),:) - coordinates(Convective(j,2),:)) * ...h(sum(coordinates(Convective(j,:),:))/2,tf,alpha)/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.figuretrisurf ( elements3, coordinates(:,1), coordinates(:,2), (full ( u )));
view(2)
colorbar
shading interp
xlabel('x')

function M = stima3 ( vertices ,K)%*****************************************************************************80
%
%% STIMA3 determines the local stiffness matrix for a triangular element.
%
%  Parameters:
%
%    Input,
%   real VERTICES(3,2), contains the 2D coordinates of the vertices.
%   K: thermal conductivity
%    Output,
%   real M(3,3), 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 );
M=M*K;

function M = stima3 ( vertices,elements ,Convective,alpha)
%*****************************************************************************80
%% STIMA3 determines the local stiffness matrix for a triangular element.
%  Parameters:
%    Input,
%   real VERTICES(3,2), contains the 2D coordinates of the vertices.
%   elements(3,1), the node index of local element
%   Convective(N,2), node index of each boundaryline of convective heat
%   transfer boundary
%   alpha, convective heat transfer coefficient
%    Output,
%   real M(3,3), the local stiffness matrix  for the element.
%
M=zeros(3,3);
for i1=1:length(Convective)temp1=ismember(elements,Convective(i1,:));if sum(temp1) == 2temp2=find(temp1 == 0);if temp2 == 1ijm=23;elseif temp2 == 2ijm=31;elseif temp2 == 3ijm=12;endswitch ijmcase 12S=norm(vertices(2,:)-vertices(1,:));M=M+alpha*S/6*[2 1 0;1 2 0;0 0 0];case 23S=norm(vertices(3,:)-vertices(2,:));M=M+alpha*S/6*[0 0 0;0 2 1;0 1 2;];case 31S=norm(vertices(3,:)-vertices(1,:));M=M+alpha*S/6*[2 0 1;0 0 0;1 0 2;];end     end
end

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) = zeros(n,1);

function value = g ( u,q )%*****************************************************************************80%% G evaluates the outward normal values assigned at Neumann boundary conditions.%  Parameters:%    Input, %    real U(N,2), contains the 2D coordinates of N points.%    q: surface heat flux at Neumann boundary%    Output, %   VALUE(N), contains the value of outward normal at each point%    where a Neumann boundary condition is applied.%value = q*ones ( size ( u, 1 ), 1 );

function value = h ( u,tf,alpha)%****************************************************************************%%  evaluates the Convective heat transfer force.%  Parameters:%    Input, %    real U(N,2), contains the 2D coordinates of N points.%    tf: ambient temperature at Convective boundary%    alpha: convective heat transfer coefficient%    Output, %   VALUE(N), contains the value of outward normal at each point%    where a Neumann boundary condition is applied.%value = alpha*tf*ones ( size ( u, 1 ), 1 );

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 = ones ( size ( u, 1 ), 1 )*100;

% 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
xminl=-10;xmaxl=10;
xminu=-5;xmaxu=5;
ymin=0;ymax=10;
Nx=13;Ny=13;
y=linspace(ymin,ymax,Ny);
xmin=linspace(xminl,xminu,Ny);
xmax=linspace(xmaxl,xmaxu,Ny);
k=0;
for i1=1:Nyx=linspace(xmin(i1),xmax(i1),Nx);for i2=1:Nxk=k+1;Coord(k,:)=[x(i2),y(i1)];    end
endsave coordinates.dat Coord -ascii% the element index
k=0;
vertices=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;vertices(k,:)=[ijm1,ijm2,ijm3];k=k+1;ijm1=i2+1+i1*Nx;ijm2=i2+i1*Nx;ijm3=i2+(i1-1)*Nx;vertices(k,:)=[ijm1,ijm2,ijm3];end
end
save elements3.dat vertices -ascii% The direchlet boundary condition (index of the two end nodes for each boundary line)
boundary=zeros(Nx-1,2);
temp1=1:Nx-1;
temp2=2:Nx;
temp3=1:Nx-1;
boundary(temp3',:)=[temp1',temp2'];
save dirichlet.dat boundary -ascii% The Neumann boundary condition (index of the two end nodes for each boundary line)
boundary=zeros(Nx-1,2);
temp1=(1:Nx-1)+(Ny-1)*Nx;
temp2=(2:Nx)+(Ny-1)*Nx;
temp3=1:Nx-1;
boundary(temp3',:)=[temp1',temp2'];
save neumann.dat boundary -ascii% The Convective heat transfer boundary condition (index of the two end nodes for each boundary line)
boundary=zeros((Ny-1)*2,2);
temp1=Nx*(1:Ny-1);
temp2=temp1+Nx;
temp3=1:Ny-1;
boundary(temp3',:)=[temp1',temp2'];
temp1=Nx*(Ny-1)+1:-Nx:Nx+1;
temp2=temp1-Nx;
temp3=temp3+Ny-1;
boundary(temp3',:)=[temp1',temp2'];
save Convective.dat boundary -ascii