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.
95 lines
2.7 KiB
95 lines
2.7 KiB
3 years ago
|
%--------------------------
|
||
|
|
% @Author: Jingqiao Hu
|
||
|
|
% @Date: 2021-01-10 11:15:50
|
||
|
|
% @LastEditTime: 2021-04-04 16:37:42
|
||
|
|
|
||
|
|
% optBodyF = 1 : Null body force
|
||
|
|
% optBodyF = 2 : Effective body force
|
||
|
|
% optBodyF = 3 : Pointwise body force
|
||
|
|
%--------------------------
|
||
|
|
function [Ue] = microFEA_QBC_2mat(D0, D1, ke0, ke1, iK, jK, edofMat, dx,dy, ...
|
||
|
|
optBodyF, ufixed, freedofs, fixeddofs, xPhys, cx, cy)
|
||
|
|
|
||
|
|
testNum = size(ufixed, 2);
|
||
|
|
e0 = eye(testNum);
|
||
|
|
F = zeros(size(ufixed));
|
||
|
|
|
||
|
|
gaussB = shapeDerivatives(dx/2, dy/2,2);
|
||
|
|
jac = dx*dy/4;
|
||
|
|
|
||
|
|
[nely, nelx] = size(xPhys);
|
||
|
|
sK = zeros(64, nelx*nely);
|
||
|
|
|
||
|
|
D = 0;
|
||
|
|
for i = 1:nelx
|
||
|
|
for j = 1:nely
|
||
|
|
ele = j + (i-1)*nely;
|
||
|
|
if xPhys(ele)==1
|
||
|
|
De = D1;
|
||
|
|
sK(:, ele) = ke1(:);
|
||
|
|
else
|
||
|
|
De = D0;
|
||
|
|
sK(:, ele) = ke0(:);
|
||
|
|
end
|
||
|
|
|
||
|
|
if optBodyF==2
|
||
|
|
D = D + De;
|
||
|
|
|
||
|
|
elseif optBodyF==3
|
||
|
|
D = De;
|
||
|
|
|
||
|
|
fe = zeros(8,testNum);
|
||
|
|
for gp = 1:4
|
||
|
|
for k = 4:testNum
|
||
|
|
eta = [e0(4, k), e0(8, k);
|
||
|
|
e0(5, k), e0(7, k);
|
||
|
|
e0(9, k), e0(6, k)] * [cx(j,i);cy(j,i)];
|
||
|
|
fe(:,k) = fe(:,k) + gaussB{gp}'*D*eta*jac;
|
||
|
|
end
|
||
|
|
end
|
||
|
|
edof = edofMat(ele,:);
|
||
|
|
F(edof,1:testNum) = F(edof,1:testNum) + fe;
|
||
|
|
end
|
||
|
|
end
|
||
|
|
end
|
||
|
|
|
||
|
|
if optBodyF==2
|
||
|
|
D = D/nelx/nely;
|
||
|
|
for i = 1:nelx
|
||
|
|
for j = 1:nely
|
||
|
|
ele = j + (i-1)*nely;
|
||
|
|
|
||
|
|
fe = zeros(8,testNum);
|
||
|
|
for gp = 1:4
|
||
|
|
for k = 4:testNum
|
||
|
|
eta = [e0(4, k), e0(8, k);
|
||
|
|
e0(5, k), e0(7, k);
|
||
|
|
e0(9, k), e0(6, k)] * [cx(j,i);cy(j,i)];
|
||
|
|
fe(:,k) = fe(:,k) + gaussB{gp}'*D*eta*jac;
|
||
|
|
end
|
||
|
|
end
|
||
|
|
|
||
|
|
edof = edofMat(ele,:);
|
||
|
|
F(edof,1:testNum) = F(edof,1:testNum) + fe;
|
||
|
|
end
|
||
|
|
end
|
||
|
|
end
|
||
|
|
|
||
|
|
K = sparse(iK, jK, sK(:));
|
||
|
|
K = (K + K') / 2;
|
||
|
|
|
||
|
|
fr = -K(freedofs, fixeddofs)*ufixed(fixeddofs,:) + F(freedofs,:);
|
||
|
|
|
||
|
|
Ue = zeros(size(ufixed));
|
||
|
|
Ue(fixeddofs, :) = ufixed(fixeddofs,:);
|
||
|
|
Ue(freedofs, :) = K(freedofs, freedofs) \ fr;
|
||
|
|
|
||
|
|
% % test
|
||
|
|
% for j = 1:testNum
|
||
|
|
% figure;
|
||
|
|
% coorx = reshape(Ue(1:2:end,j),nely+1,nelx+1) + cx;
|
||
|
|
% coory = reshape(Ue(2:2:end,j),nely+1,nelx+1) + cy;
|
||
|
|
% scatter(coorx(:),coory(:)); hold on; axis equal;
|
||
|
|
% end
|
||
|
|
end
|