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.
85 lines
2.6 KiB
85 lines
2.6 KiB
4 years ago
|
% optBodyF = 1 : Null body force
|
||
|
|
% optBodyF = 2 : Effective body force
|
||
|
|
% optBodyF = 3 : Pointwise body force
|
||
|
|
function [Ue] = microFEA_QBC(De, ke, iK, jK, edofMat, dx,dy, ...
|
||
|
|
optBodyF, ufixed, freedofs, fixeddofs, xPhys, cx,cy, e, emin)
|
||
|
|
|
||
|
|
testNum = size(ufixed, 2);
|
||
|
|
e0 = eye(testNum);
|
||
|
|
F = zeros(size(ufixed));
|
||
|
|
gaussB = shapeDerivatives(dx,dy);
|
||
|
|
jac = dx*dy/4;
|
||
|
|
|
||
|
|
[nely, nelx] = size(xPhys);
|
||
|
|
|
||
|
|
D = 0;
|
||
|
|
for i = 1:nelx
|
||
|
|
for j = 1:nely
|
||
|
|
ele = j + (i-1)*nely;
|
||
|
|
De = De * xPhys(j,i);
|
||
|
|
|
||
|
|
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
|
||
|
|
|
||
|
|
if nargin>14
|
||
|
|
sK = ke(:) * (emin + xPhys(:)'*(e-emin));
|
||
|
|
else
|
||
|
|
sK = ke(:) * xPhys(:)'; % xPhys = sum(H(De)^2)/4, thus not need penalty
|
||
|
|
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;
|
||
|
|
% end
|
||
|
|
end
|