From d20381238b897747946f09f0c03e0ed77a7641f3 Mon Sep 17 00:00:00 2001
From: "lisa.pizzo" <lisa.pizzo@noreply.imsc.uni-graz.at>
Date: Tue, 27 Jan 2026 09:47:35 +0100
Subject: [PATCH] Task 4

---
 Project/ApplyRobinBC_mult_rot.m | 36 +++++++++++++++++++++++++++++++++
 1 file changed, 36 insertions(+)
 create mode 100644 Project/ApplyRobinBC_mult_rot.m

diff --git a/Project/ApplyRobinBC_mult_rot.m b/Project/ApplyRobinBC_mult_rot.m
new file mode 100644
index 0000000..3a9d8ba
--- /dev/null
+++ b/Project/ApplyRobinBC_mult_rot.m
@@ -0,0 +1,36 @@
+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);
+    
+    %find the edges on the axis -> homogeneous Neumann BC, no Robin contribution
+    if ri == 0 && rj == 0
+        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