% Solution to Challenge 2, Homework 4, CMSC/AMSC 661,
% Spring 2010  Dianne P. O'Leary
% This problem is inspired by one by Gunilla Kreiss, lab6.ps
% http://www.nada.kth.se//kurser/kth/2D1250/tilnum2-04/index.html

global k_copper k_silicon q_silicon h_heat_transf

k_copper  = .400;     % units: W/(mm K)
k_silicon = .148;     % units: W/(mm K)
q_silicon = 5;  % units: W/(mm^3)  (original: 5.86e7)
h_heat_transf = 1.e-5; % units: W/(mm^2 K)
room_temp = 300;       % units: K

% Step 1: set up the x-y mesh.
%
% The mesh is specified by the %  --  Mesh points, stored in p.
%       (A mesh with n points makes p 2xn.  The x-coordinates
%        are in row 1 and the y-coordinates are in row 2.)
%  --  Triangle edge information, stored in e.
%       (This array is not so easy to decode.)
%  --  Triangle vertex information, stored in t.
%       (If column 10 of t is [212, 43, 5, 1]', then
%        the 10th triangle has mesh points 212, 43, and 5
%        as its vertices.)

% We define a rectangular domain with vertices (0,0) and (20,4+a).
% The matrix d1 can be determined by the Matlab statements
%   pderect([0 20 0 4])
%   < Then click on Draw -> Export  >
%   [d1,bt]=decsg(gd,sf,ns);
% The first column indicates that the boundary segment is a
% straight line (d1(1,1)=1) with x-range 0 to 20 and y-range
% 0 to 0.

jj = 0;
thickness = [linspace(0,10,11),linspace(4.0,4.2,5)];

for myb = thickness
hght = 2 + myb;
d1 = [ 2     2     2     2
       0    20    20     0
      20    20     0     0
       0     0  hght  hght
       0  hght  hght     0
       1     1     1     1
       0     0     0     0 ];

bt = 1;

% In this example, we let Matlab divide the region into a
% coarse mesh of triangles, and then refine the mesh once.
% It would be much better to specify the domain so that
% no triangle overlaps the silicon and the copper regions.

[p,e,t] = initmesh(d1);
[p,e,t] = refinemesh(d1,p,e,t);
[p,e,t] = refinemesh(d1,p,e,t);
[p,e,t] = refinemesh(d1,p,e,t);
[p,e,t] = refinemesh(d1,p,e,t);
% figure(1)
% pdemesh(p,e,t)

% Step 2: Specify initial conditions.
%
% Now we set the initial conditions v(x1,x2) for our pde. 
% In this example, we set the initial temperature to be
% room temperature.
%

u0 = room_temp*ones(size(p,2),1);

% Step 3: Set the times of interest.
%
% Now we specify the times at which we want to compute
% the solution: 0, 1, ..., 100.              

tlist = linspace(0, 100, 101);

% Step 4a: Specify the boundary conditions.
%
% This is not well documented by Mathworks, but a good
% example can be found in:
% http://www.mathworks.com/matlabcentral/fx_files/9824/1/solveBV.m
%
% Our boundary conditions are du/dn = 0 for the bottom
% boundary, and -c du/dn = h (u-u*) for the other boundaries,
% where u* = room temperature % This is specified by the matrix b.
% Matlab assumes that the boundary conditions are in the
% format  n . (c grad u) + q u = g,
% so for our problem, q = h and g = h*room_temp.

z  = mat2str(0);
tmp = h_heat_transf*room_temp;
g = mat2str(tmp);
q  = mat2str(h_heat_transf);
lz = length(z);
lg = length(g);
lq = length(q);
z1q = z*ones(1,lq);
z1g = z*ones(1,lg);
bc = [1      1       1       1
      0      0       0       0
      lz     lq      lq      lq
      lz     lg      lg      lg
       z1q'  (q*1)' (q*1)' (q*1)'
       z1g'  (g*1)' (g*1)' (g*1)' ];

% Step 4b. Specify the pde coefficients.
%
% The last four arguments to Matlab's parabolic
%     are c, a, f, and d.  These specify the pde to be
%      d u_t - grad . (c grad u) + a u = f.

d = 1;
a = 0;

% We need to set the coefficients c and f at the centers of
% mass of each triangle.
% Define copper(j) = 1 if the center of mass of triangle j
% is in the copper region, and 0 otherwise.

[kk,ntri] = size(t);

is_copper = [];

for j=1:ntri,
   coord = (p(:,t(1,j))+p(:,t(2,j))+p(:,t(3,j)) ) / 3;
   if ( (coord(1) < 2) | (coord(1) > 18) | (coord(2) > 2) )
      is_copper(j) = 1;
   else
      is_copper(j) = 0;
   end
end

c = (1-is_copper)*k_silicon + is_copper*k_copper;
f = (1-is_copper)*q_silicon;

% Step 4c. Solve the pde.

u1 = parabolic(u0, tlist, bc, p, e, t, c, a, f, d);

% Step 5: Use the solution.
%
% u1(i,j) is the approximate solution at the i-th
% mesh point at the j-th time.
% The next commands plot the solution at time 0.01
% and at the final time.

%figure(2)
%pdeplot(p,e,t,'xydata',u1(:,2))

%figure(3)
%pdeplot(p,e,t,'xydata',u1(:,end))

jj = jj + 1;
maxu(jj) = max(u1(:,end));
end
[thickness;maxu]'
pdeplot(p,e,t,'xydata',u1(:,end))
