lisa.pizzo/LisaPizzoExercises
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Update Project/ApplyRobinBC_mult_rot.m

Browse source
This commit is contained in:
Lisa Pizzo 2026-01-28 13:13:58 +01:00
commit d69f21fb6f
1 changed files with 29 additions and 28 deletions
Download patch file Download diff file Expand all files Collapse all files

35
Project/ApplyRobinBC_mult_rot.m
View file

@ -1,36 +1,37 @@
function [K, F] = ApplyRobinBC_mult_rot(model, K, F, alpha, u_out) function [K, F] = ApplyRobinBC_mult_rot(model, K, F, alpha, u_out)
mesh = model.Mesh; mesh = model.Mesh;
nodes = mesh.Nodes; nodes = mesh.Nodes;
elements = mesh.Elements;
edgesAll = [elements([1 2],:), elements([2 3],:), elements([3 1],:)]; robinEdges = [4, 5, 6, 11];
edgesSorted = sort(edgesAll,1);
[edgesUnique,~,ic] = unique(edgesSorted','rows');
counts = accumarray(ic,1);
boundaryEdges = edgesUnique(counts==1,:);
for k = 1:size(boundaryEdges,1) for eID = robinEdges
i = boundaryEdges(k,1);
j = boundaryEdges(k,2);
ri = nodes(1,i); edgeNodes = findNodes(mesh,'region','Edge',eID);
rj = nodes(1,j); for k = 1:length(edgeNodes)-1
i = edgeNodes(k);
%find the edges on the axis -> homogeneous Neumann BC, no Robin contribution j = edgeNodes(k+1);
if ri == 0 && rj == 0
continue;
end
xi = nodes(:,i); xi = nodes(:,i);
xj = nodes(:,j); xj = nodes(:,j);
% Edge length
L = norm(xi - xj); L = norm(xi - xj);
rbar = 0.5 * (xi(1) + xj(1)); % r at edge midpoint
% Axisymmetric radius at midpoint
rbar = 0.5 * (xi(1) + xj(1));
% Skip axis (safety, should not occur here)
if rbar == 0
continue;
end
% Robin element matrices (linear edge)
Ke = rbar * alpha * L / 6 * [2 1; 1 2]; Ke = rbar * alpha * L / 6 * [2 1; 1 2];
Fe = rbar * alpha * u_out * L / 2 * [1; 1]; Fe = rbar * alpha * u_out * L / 2 * [1; 1];
% Assemble
K([i j],[i j]) = K([i j],[i j]) + Ke; K([i j],[i j]) = K([i j],[i j]) + Ke;
F([i j]) = F([i j]) + Fe; F([i j]) = F([i j]) + Fe;
end
end end
end end