a 2D version
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 

126 lines
3.9 KiB

% Quadratic boundary conditions
function [ufixed,border, freedofs] = microQBC(nelx,nely,cx,cy)
nodenrs = reshape(1:(1 + nelx) * (1 + nely), 1 + nely, 1 + nelx);
e0 = eye(9);
testNum = size(e0,1);
ufixed = zeros(2*(nely + 1)*(nelx + 1), testNum);
alldofs = 1:2 * (nely + 1) * (nelx + 1);
n1 = [nodenrs(end, [1, end]), nodenrs(1, [end, 1])];
d1 = reshape([2 * n1 - 1; 2 * n1], 1, 8); % four corner points
n3 = [nodenrs(2:end - 1, 1)', nodenrs(end, 2:end - 1)];
d3 = reshape([2 * n3 - 1; 2 * n3], 1, 2 * (nelx + nely - 2)); % left+bottom boundaries
n4 = [nodenrs(2:end - 1, end)', nodenrs(1, 2:end - 1)];
d4 = reshape([2 * n4 - 1; 2 * n4], 1, 2 * (nelx + nely - 2)); % right+up boundaries
d2 = setdiff(alldofs, [d1, d4, d3]);
coor1m = [cx(end, [1, end]), cx(1, [end, 1]);
cy(end, [1, end]), cy(1, [end, 1])];
coor3m = [cx(2:end - 1, 1)', cx(end, 2:end - 1);
cy(2:end - 1, 1)', cy(end, 2:end - 1)];
coor4m = [cx(2:end - 1, end)', cx(1, 2:end - 1);
cy(2:end - 1, end)', cy(1, 2:end - 1)];
coor1 = coor1m(:);
coor3 = coor3m(:);
coor4 = coor4m(:);
for j = 1:3
e = [e0(1, j), e0(3, j) / 2; e0(3, j) / 2, e0(2, j)];
% d1, corner
for i = 1:2:length(d1)
ufixed(d1(i:i+1), j) = e * coor1(i:i+1,1);
end
% border
for i = 1:2:length(d3)
ufixed(d3(i:i+1), j) = e * coor3(i:i+1,1);
ufixed(d4(i:i+1), j) = e * coor4(i:i+1,1);
end
end
if testNum>3
for j = 4:testNum
[ex,ey] = strainGradient(e0(:,j));
% d1, corner
for i = 1:2:length(d1)
coor = [coor1(i,1), coor1(i+1,1)];
ufixed(d1(i), j) = 1/2 * coor * ex * coor';
ufixed(d1(i+1), j) = 1/2 * coor * ey * coor';
end
% d3, left+bottom
for i = 1:2:length(d3)
coor = [coor3(i,1),coor3(i+1,1)];
ufixed(d3(i), j) = 1/2 * coor * ex * coor';
ufixed(d3(i+1), j) = 1/2 * coor * ey * coor';
coor = [coor4(i,1),coor4(i+1,1)];
ufixed(d4(i), j) = 1/2 * coor * ex * coor';
ufixed(d4(i+1), j) = 1/2 * coor * ey * coor';
end
end
end
border = [d1,d3,d4];
freedofs = d2;
% test
% for j = 1:testNum
% figure;
% for i = 1:2:length(d1)
% coor = ufixed(d1(i:i+1),j)+coor1(i:i+1,1);
% scatter(coor(1:2:end),coor(2:2:end),'filled'); hold on;
% end
%
% for i = 1:2:length(d3)
% coor = ufixed(d3(i:i+1), j) + coor3(i:i+1,1);
% scatter(coor(1:2:end),coor(2:2:end)); hold on;
% coor = ufixed(d4(i:i+1), j) + coor4(i:i+1,1);
% scatter(coor(1:2:end),coor(2:2:end)); hold on;
% end
% axis equal;
% end
end
function u = bilinearU(x,y,lx,ly,u0)
N1 = (1-x/lx).*(1-y/ly);
N2 = x/lx.*(1-y/ly);
N3 = x/lx.*y/ly;
N4 = (1-x/lx).*y/ly;
u1 = N1*u0(1) + N2*u0(3) + N3*u0(5) + N4*u0(7);
u2 = N1*u0(2) + N2*u0(4) + N3*u0(6) + N4*u0(8);
u = [u1,u2];
end
function [gx,gy] = strainGradient(e0)
e = zeros(2,2,2);
e(1,1,1) = e0(4);
e(2,2,1) = e0(5);
e(1,2,2) = e0(6); e(2,1,2) = e(1,2,2);
e(2,2,2) = e0(7);
e(1,1,2) = e0(8);
e(1,2,1) = e0(9); e(2,1,1) = e(1,2,1);
g111 = secondGradient(e,1,1,1);
g121 = secondGradient(e,1,2,1);
g112 = secondGradient(e,1,1,2);
g122 = secondGradient(e,1,2,2);
gx = [g111, g121;
g112, g122];
g211 = secondGradient(e,2,1,1);
g221 = secondGradient(e,2,2,1);
g212 = secondGradient(e,2,1,2);
g222 = secondGradient(e,2,2,2);
gy = [g211, g221;
g212, g222];
end
function g = secondGradient(e,i,j,k)
g = e(i,j,k)+e(i,k,j) - e(j,k,i);
end