%% 2d heat conduction
%
% Numerical Analysis of Structures WS 2017
%
% (c) 2017 Marek Tyburec, marek.tyburec@fsv.cvut.cz

%% DEFINE STRUCTURE
% coordinates of the nodes
xCoordinates = [0 0.3 0.3 0.3 0.38 0.38 0 0.58 0.58 0 0.0375 0.075 0.1125 0.15 0.1875 0.225 0.2625 0.3 0.3 0.3 0.2625 0.225 0.1875 0.15 0.1125 0.075 0.0375 0 0 0 0.3 0.34 0.38 0.342 0.304 0.266 0.228 0.19 0.152 0.114 0.076 0.038 0 0 0 0.34 0.38 0.42 0.46 0.5 0.54 0.58 0.54 0.5 0.46 0.42 0.023978596 0.053670673 0.091222877 0.13055405 0.16948683 0.20878494 0.24630498 0.27598145 0.27130555 0.27131845 0.27598005 0.24629582 0.20876546 0.16942855 0.13049622 0.091196022 0.053679683 0.024011604 0.028688216 0.028722306 0.065459974 0.10992766 0.14998654 0.19003981 0.23449821 0.27591002 0.27076329 0.32130152 0.35721508 0.28459407 0.24859154 0.21055938 0.1710594 0.13183355 0.09216709 0.054274576 0.024287125 0.028945224 0.028845188 0.024089907 0.053877471 0.091546936 0.13100024 0.16996058 0.20909294 0.2462868 0.26768239 0.23478001 0.19104342 0.15099891 0.11074042 0.065978697 0.55304171 0.54546441 0.32695834 0.33453566 0.36025602 0.40004266 0.44000001 0.47995737 0.51974413];
yCoordinates = [0 0 0.15 0.23 0.23 0.3 0.3 0.15 0.23 0.15 0 0 0 0 0 0 0 0.0375 0.075 0.1125 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.1125 0.075 0.0375 0.19 0.23 0.265 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.2625 0.225 0.1875 0.15 0.15 0.15 0.15 0.15 0.15 0.19 0.23 0.23 0.23 0.23 0.025149945 0.032033028 0.036320151 0.037000379 0.037012753 0.036346782 0.032057292 0.025156437 0.056297567 0.09346649 0.12470958 0.11766142 0.11320098 0.11239532 0.11236602 0.11312674 0.11759324 0.12466853 0.093396293 0.056216697 0.074594285 0.074464035 0.074449694 0.074537005 0.074696979 0.17586424 0.20691563 0.26527415 0.26505037 0.26659454 0.26546837 0.26300008 0.26280879 0.26293372 0.26363871 0.26793904 0.27484259 0.24370007 0.20653099 0.17528932 0.18233427 0.1867834 0.1875563 0.18750348 0.18660771 0.18240559 0.23879666 0.22426703 0.22516288 0.22543334 0.22542828 0.22529456 0.17604156 0.20660814 0.20395856 0.17339191 0.19640086 0.19106646 0.18999975 0.18893322 0.18359968];
numberOfNodes = length(xCoordinates);

% elements (node numbers counterclockwise)
elements = [1 11 57; 11 12 58; 12 13 59; 13 14 60; 14 15 61; 15 16 62; 16 17 63; 17 2 64; 2 18 64; 18 19 65; 19 20 66; 20 3 67; 3 21 67; 21 22 68; 22 23 69; 23 24 70; 24 25 71; 25 26 72; 26 27 73; 27 10 74; 10 28 74; 28 29 75; 29 30 76; 30 1 57; 57 11 58; 58 12 59; 59 13 60; 60 14 61; 61 15 62; 62 16 63; 63 17 64; 64 18 65; 65 19 66; 66 20 67; 67 21 68; 68 22 69; 69 23 70; 70 24 71; 71 25 72; 72 26 73; 73 27 74; 74 28 75; 75 29 76; 76 30 57; 57 58 76; 58 59 77; 59 60 78; 60 61 79; 61 62 80; 62 63 81; 63 64 65; 65 66 81; 66 67 68; 68 69 81; 69 70 80; 70 71 79; 71 72 78; 72 73 77; 73 74 75; 75 76 77; 76 58 77; 77 59 78; 78 60 79; 79 61 80; 80 62 81; 81 63 65; 81 66 68; 81 69 80; 80 70 79; 79 71 78; 78 72 77; 77 73 75; 3 31 82; 31 4 83; 4 32 84; 32 5 85; 5 33 85; 33 6 85; 6 34 85; 34 35 84; 35 36 86; 36 37 87; 37 38 88; 38 39 89; 39 40 90; 40 41 91; 41 42 92; 42 7 93; 7 43 93; 43 44 94; 44 45 95; 45 10 96; 10 27 96; 27 26 97; 26 25 98; 25 24 99; 24 23 100; 23 22 101; 22 21 102; 21 3 82; 82 31 83; 83 4 103; 4 84 86; 84 32 85; 85 34 84; 84 35 86; 86 36 87; 87 37 88; 88 38 89; 89 39 90; 90 40 91; 91 41 92; 92 42 93; 93 43 94; 94 44 95; 95 45 96; 96 27 97; 97 26 98; 98 25 99; 99 24 100; 100 23 101; 101 22 102; 102 21 82; 82 83 102; 83 103 104; 103 4 86; 86 87 103; 87 88 104; 88 89 105; 89 90 106; 90 91 107; 91 92 108; 92 93 94; 94 95 108; 95 96 97; 97 98 108; 98 99 107; 99 100 106; 100 101 105; 101 102 104; 102 83 104; 104 103 87; 104 88 105; 105 89 106; 106 90 107; 107 91 108; 108 92 94; 108 95 97; 108 98 107; 107 99 106; 106 100 105; 105 101 104; 8 52 109; 52 9 110; 4 31 111; 31 3 112; 5 32 113; 32 4 111; 3 46 112; 46 47 112; 47 48 114; 48 49 115; 49 50 116; 50 51 117; 51 8 109; 9 53 110; 53 54 110; 54 55 116; 55 56 115; 56 5 114; 109 52 110; 111 31 112; 5 113 114; 113 32 111; 112 47 113; 47 114 113; 114 48 115; 115 49 116; 116 50 117; 117 51 109; 110 54 117; 54 116 117; 116 55 115; 115 56 114; 109 110 117; 111 112 113];
numberOfElements = size(elements,1);

% coefficient of thermal conductivity
lambda = zeros(numberOfElements, 1);
lambda(1:72) = 0.33;
lambda(73:152) = 0.038;
lambda(153:end) = 0.104;

% inner source of heat (constant over element)
sourceElements = [];
sourceValue = [];

% heat flux
fluxElements = [];
fluxEdge = [];
fluxValue = [];

% heat transfer
transferElements = [1 2 3 4 5 6 7 8 9 10 11 12 159 160 161 162 163 164 165 166 167 168 169 170 77 78 79 80 81 82 83 84 85 86 87 88];
transferEdge = ones(1, length(transferElements));
transferTemperature = zeros(1, length(transferElements));
transferTemperature(1:19) = 20;
transferTemperature(20:end) = -10;
transferFilmCoefficient = zeros(1, length(transferElements));
transferFilmCoefficient(1:19) = 7.7;
transferFilmCoefficient(20:end) = 25;

% prescribed temperature
nodeNumberTemperature = [];
temperature = [];

%% FEM ANALYSIS
% preallocation
conductivityMatrix = zeros(numberOfNodes, numberOfNodes);
temperatureVector = zeros(numberOfNodes, 1);
rhs = zeros(numberOfNodes,1);
areas = zeros(numberOfElements,1);

% conductivity matrix assembly
for element = 1:numberOfElements
    % pick the element 'element'
    
    % nodal coordinates
    x1 = xCoordinates(elements(element, 1));
    x2 = xCoordinates(elements(element, 2));
    x3 = xCoordinates(elements(element, 3));
    y1 = yCoordinates(elements(element, 1));
    y2 = yCoordinates(elements(element, 2));
    y3 = yCoordinates(elements(element, 3));
    
    % area
    directionVector1 = [x2-x1, y2-y1];
    directionVector2 = [x3-x1, y3-y1];
    areas(element) = (directionVector1(1)*directionVector2(2) - directionVector1(2)*directionVector2(1))/2;
    
    % ----------------------------------------------------------------------------------------------------
    % CREATE ELEMENT CONDUCTIVITY MATRIX (conductivityElement variable)
    % ----------------------------------------------------------------------------------------------------
    
    % code numbers
    codeNumbers = elements(element,:);
    
    % assembly
    conductivityMatrix(codeNumbers, codeNumbers) = ...
        conductivityMatrix(codeNumbers, codeNumbers) + conductivityElement;
end

% right hand side for inner source of heat
for source = 1:length(sourceElements)
    element = sourceElements(source);
    codeNumbers = elements(element,:);
    rhs(codeNumbers) = rhs(codeNumbers) + sourceValue(source) * areas(element) / 3 * [1; 1; 1];
end

% right hand side for prescribed heat flux
nodeOrder = [1 2 3 1];
for flux = 1:length(fluxElements)
    element = fluxElements(flux);
    codeNumbers = elements(element,:);
    firstNodeOrder = nodeOrder(fluxEdge(flux));
    secondNodeOrder = nodeOrder(fluxEdge(flux)+1);
    edgeDX = xCoordinates(elements(element, secondNodeOrder) ) - ...
             xCoordinates(elements(element, firstNodeOrder) );
    edgeDY = yCoordinates(elements(element, secondNodeOrder) ) - ...
             yCoordinates(elements(element, firstNodeOrder ) );
    edgeLength = sqrt(edgeDX.^2 + edgeDY.^2);
    nodalDistribution = zeros(3,1);
    nodalDistribution([firstNodeOrder secondNodeOrder]) = 1;
    rhs(codeNumbers) = rhs(codeNumbers) - fluxValue(flux) * edgeLength/2 * nodalDistribution;
end

% heat transfer (right hand side + convective stiffness matrix)
nodeOrder = [1 2 3 1];
for transfer = 1:length(transferElements)
    element = transferElements(transfer);
    codeNumbers = elements(element,:);
    firstNodeOrder = nodeOrder(transferEdge(transfer));
    secondNodeOrder = nodeOrder(transferEdge(transfer)+1);
    edgeDX = xCoordinates(elements(element, secondNodeOrder) ) - ...
             xCoordinates(elements(element, firstNodeOrder) );
    edgeDY = yCoordinates(elements(element, secondNodeOrder) ) - ...
             yCoordinates(elements(element, firstNodeOrder ) );
    edgeLength = sqrt(edgeDX.^2 + edgeDY.^2);
    nodalDistribution = zeros(3,1);
    nodalDistribution([firstNodeOrder secondNodeOrder]) = 1;
    
    % element convective matrix preallocation
    conductivityElement = zeros(3,3);
    conductivityElement(find(nodalDistribution), find(nodalDistribution)) = ...
        transferFilmCoefficient(transfer) * edgeLength * [1/3 1/6; 1/6 1/3];
    
    % assembly
    conductivityMatrix(codeNumbers, codeNumbers) = ...
        conductivityMatrix(codeNumbers, codeNumbers) + conductivityElement;
    
    % right hand side
    rhs(codeNumbers) = rhs(codeNumbers) + ...
        transferTemperature(transfer) * transferFilmCoefficient(transfer) * ...
        edgeLength / 2 * nodalDistribution;
end

% right hand side from prescribed temperature
temperatureVector(nodeNumberTemperature) = temperature;
boundaryConditionsEssential = nodeNumberTemperature;
boundaryConditionsFree = setdiff(1:numberOfNodes, boundaryConditionsEssential);
if isempty(boundaryConditionsEssential)
    rhsTemperatureEffect = zeros(length(boundaryConditionsFree),1);
else
    rhsTemperatureEffect = conductivityMatrix(boundaryConditionsFree,boundaryConditionsEssential)* ...
                                          temperatureVector(boundaryConditionsEssential);
end

% solution
temperatureVector(boundaryConditionsFree) = conductivityMatrix(boundaryConditionsFree, boundaryConditionsFree) ...
                                            \ (rhs(boundaryConditionsFree) - rhsTemperatureEffect);

%% PLOTTING
axis equal;
xlabel('x [m]'); ylabel('y [m]'); colorbar;
hold on;

% create outline
outlineX = [0 0.3 0.3 0.58 0.58 0.38 0.38 0];
outlineY = [0 0 0.15 0.15 0.23 0.23 0.3 0.3];

% resample points
[Xr, Yr, Wr] = griddata(xCoordinates, yCoordinates, temperatureVector, ...
    unique(xCoordinates),unique(yCoordinates)');

% remove outer
% note: octave plots nonsense with nans
IN = inpolygon(Xr, Yr, outlineX, outlineY);
Wr(~IN) = 0;
Wr(isnan(Wr)) = 0;

% plot
contourf(Xr, Yr, Wr, 30);

% hide triangles in non-convex regions
patch([0.3 0.58 0.58 0.3], [0 0 0.15 0.15], 1, 'FaceColor', [1 1 1], 'EdgeColor', 'none');
patch([0.38 0.58 0.58 0.38], [0.23 0.23 0.3 0.3], 1, 'FaceColor', [1 1 1], 'EdgeColor', 'none');

% overlay elements
% for element = 1:numberOfElements    
%     xElement = xCoordinates(elements(element,:));
%     yElement = yCoordinates(elements(element,:));
%     tElement = temperatureVector(elements(element,:));
%     patch(xElement([1 2 3 1]), yElement([1 2 3 1]), tElement([1 2 3 1]), 'FaceColor', 'none', 'LineStyle', '--', 'EdgeColor', 0.5*[1 1 1]);
% end

