Task 5

Project/ApplyRobinBC_mult_rot.m

@@ -0,0 +1,37 @@
+function [K, F] = ApplyRobinBC_mult_rot(model, K, F, alpha, u_out)
+mesh = model.Mesh;
+nodes = mesh.Nodes;
+elements = mesh.Elements;
+
+edgesAll = [elements([1 2],:), elements([2 3],:), elements([3 1],:)];
+edgesSorted = sort(edgesAll,1);
+[edgesUnique,~,ic] = unique(edgesSorted','rows');
+counts = accumarray(ic,1);
+boundaryEdges = edgesUnique(counts==1,:);
+
+for k = 1:size(boundaryEdges,1)
+    i = boundaryEdges(k,1);
+    j = boundaryEdges(k,2);
+
+    ri = nodes(1,i);
+    rj = nodes(1,j);
+    
+    if ri == 0 && rj == 0
+    % Edge lies on symmetry axis r = 0
+    % -> homogeneous Neumann BC, no Robin contribution
+        continue;
+    end
+
+    xi = nodes(:,i);
+    xj = nodes(:,j);
+
+    L = norm(xi - xj);
+    rbar = 0.5 * (xi(1) + xj(1));   % r at edge midpoint
+
+    Ke = rbar * alpha * L / 6 * [2 1; 1 2];
+    Fe = rbar * alpha * u_out * L / 2 * [1; 1];
+
+    K([i j],[i j]) = K([i j],[i j]) + Ke;
+    F([i j]) = F([i j]) + Fe;
+end
+end

Project/CalculateLaplace_mult_rot.m

@@ -0,0 +1,40 @@
+function [K, F] = CalculateLaplace_mult_rot(model, lambda_wall, lambda_fluid, lambda_air)
+mesh = model.Mesh;
+nodes = mesh.Nodes;
+elements = mesh.Elements;
+
+Nnodes = size(nodes,2);
+Nelems = size(elements,2);
+
+K = sparse(Nnodes, Nnodes);
+F = zeros(Nnodes,1);
+
+regions = zeros(Nelems,1);
+regions(findElements(mesh,'region','Face',1)) = 1;
+regions(findElements(mesh,'region','Face',2)) = 2;
+regions(findElements(mesh,'region','Face',3)) = 3;
+
+for e = 1:Nelems
+    vert = elements(:,e);
+
+    x = nodes(1,vert);   % r-coordinates
+    y = nodes(2,vert);   % z-coordinates
+
+    Ae = polyarea(x,y);
+
+    b = [y(2)-y(3); y(3)-y(1); y(1)-y(2)];
+    c = [x(3)-x(2); x(1)-x(3); x(2)-x(1)];
+
+    rbar = mean(x);      % <-- axisymmetric weight
+    switch regions(e)
+        case 1
+            lambda = lambda_wall;
+        case 2
+            lambda = lambda_fluid;
+        case 3
+            lambda = lambda_air;
+    end
+    Ke = rbar * (lambda/(4*Ae)) * (b*b.' + c*c.');
+    K(vert,vert) = K(vert,vert) + Ke;
+end
+end