clear all; clc
%Ref: HTn.17, "Fun" Problem, Heat Conduction in Micro-Electronics

%set physical parameters
k_a = 1;
k_b =  .1;
s_b = 200;
T_b = 100;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set up the mesh and element connectivity tables
%  -- generate the mesh points for the elipse
%  -- generate the mesh points for the circluar sink area
%  -- generate the triangular elements from the mesh points
%  -- segregate the elements into those belonging to omega_a and those
%        belonging to omega_b

nrad = 5; % nodes in radius direction
ntheta = 20; % nodes in theta direction
a=1;
b=.75;
rsweep1=linspace(.1,1,nrad);
rsweep2=linspace(.1,1,2*nrad);
dtheta = 2*pi/ntheta;
X1= zeros(1,ntheta*4*nrad);   % reserve memory
Y1= zeros(1,ntheta*4*nrad);   % reserve memory

theta = 0;
ptr =0;
for i=1:ntheta
  r = 1/sqrt((cos(theta)/a)^2 + (sin(theta)/b)^2);
  if theta>pi/2 & theta<3*pi/2;
    rsweep = rsweep1;
  else 
    rsweep = rsweep2;
  end
  rads = r * rsweep;
  nrads = length(rads);
  X1(ptr+1:ptr+nrads) = rads * cos(theta);
  Y1(ptr+1:ptr+nrads) = rads * sin(theta);
  ptr = ptr + nrads ;
  theta = theta + dtheta;
end

% generate extra points near .422, -.222
r=.2;
rads = r*rsweep2;
nrads = length(rads);
theta = 0;
for i=1:ntheta
  X1(ptr+1:ptr+nrads) = .422 + rads * cos(theta);
  Y1(ptr+1:ptr+nrads) = -.222 +  rads * sin(theta);
  ptr = ptr + nrads ;
  theta = theta + dtheta;
end

% generate a point at .422, -.222
X1(ptr+1) = .422;
Y1(ptr+1) = -.222;
ptr = ptr + 1;

% trim zeros off the end
X=X1(1:ptr+1);
Y=Y1(1:ptr+1);

% use Matlab's implementation of the Delaunay trianglation algorithm to
% generate the triangles and the corresponding element connectivity table
tri = delaunay(X,Y);
elcon = fliplr(tri); % reverse the order of nodes in each element

% uncomment to view the mesh
 trimesh(elcon,X,Y,zeros(size(X))); view(2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculate the center point of each triangle, then sort the elements
% based on whether or not the center lies inside or outside of the
% circular domain

X_center = sum(X(elcon),2)/3;
Y_center = sum(Y(elcon),2)/3;
% calculate the distance from the center of each element to the center of 
% omega2
distance = sqrt((X_center - .422).^2 + (Y_center + .222).^2);

% initialize element conductivity tables for omega1 and omega2
elcon_omega_a = zeros(size(elcon));
elcon_omega_b = zeros(size(elcon));
ptr_a = 1; % pointer to next free position in elcon_omega_a
ptr_b = 1; % pointer to next free position in elcon_omega_b

% cycle through the elements (rows of elcon) and copy rows to
% elcon_omega_a or elcon_omega_b.
for i=1:size(elcon,1)
  if distance(i) > .2
    elcon_omega_a(ptr_a,:) = elcon(i,:);
    ptr_a = ptr_a + 1;
  else
    elcon_omega_b(ptr_b,:) = elcon(i,:);
    ptr_b = ptr_b + 1;
  end
end
% trim off unused rows from elcon_omega_a and elcon_omega_b
elcon_omega_a = elcon_omega_a(1:ptr_a-1,:);
elcon_omega_b = elcon_omega_b(1:ptr_b-1,:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% mesh and element connectivity tables are done, now assemble the
% elements

% cacluate ZETAe, DETe, and DIFFe for all the elements in omega_a and omega_b
[ZETAe_omega_a,DETe_omega_a]=do_ZETAe_2Dt(X,Y,elcon_omega_a);
[ZETAe_omega_b,DETe_omega_b]=do_ZETAe_2Dt(X,Y,elcon_omega_b);
DIFFe_omega_a = do_DIFFe_2Dt(ZETAe_omega_a); % element diffusion matrices
DIFFe_omega_b = do_DIFFe_2Dt(ZETAe_omega_b); % element diffusion matrices

% assemble diffusion matrices 
DIFF_omega_a =  asjac(k_a,[],[],-1,DIFFe_omega_a,[],...
		      elcon_omega_a,DETe_omega_a,prod(size(X)));
DIFF_omega_b =  asjac(k_b,[],[],-1,DIFFe_omega_b,[],...
		      elcon_omega_b,DETe_omega_b,prod(size(X)));

% assemble the source term for omega_b.  
B200L = [2 1 1; ...
	 1 2 1; ...
	 1 1 2]/24;
SRC = s_b * ones(size(X(:)));
b = asres([],[],[],1,B200L,SRC,...
	  elcon_omega_b,DETe_omega_b,prod(size(X)));

DIFF_omega = DIFF_omega_a + DIFF_omega_b;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% apply dirichlet BC to all points on the outer edge of omega.
% first, make a table of the nodes to have diri BC applied
diri_nodes = convhull(X,Y);
diri_list = zeros(prod(size(X)),2);
for i=1:prod(size(diri_nodes))
  diri_list(diri_nodes(i),:) = [1 T_b];
end
[DIFF_diri, b_diri] = do_diri(DIFF_omega,b,diri_list);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve the linear system
Q = DIFF_diri \  b_diri;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot the solution
trimesh(tri,X,Y,Q); view(3);