hjq
/
Numerical_coarsening_using_discontinuous_shape_functions
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
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.
2 Commits
1 Branch
0 Tags
707 KiB
Branch: master
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from 'master'
${ noResults }
Numerical_coarsening_using_.../utility-func/getStress.m

126 lines
14 KiB
Raw Permalink Normal View History

init
3 years ago
clear all
% lx = 1; ly = 1;
syms x y lx ly
N(1) = 1/4*(1-x/lx)*(1-y/ly);
N(2) = 1/4*(1+x/lx)*(1-y/ly);
N(3) = 1/4*(1+x/lx)*(1+y/ly);
N(4) = 1/4*(1+x/lx)*(1+y/ly);
Be= [diff(N(1),x) 0 diff(N(2),x) 0 diff(N(3),x) 0 diff(N(4),x) 0;
0 diff(N(1),y) 0 diff(N(2),y) 0 diff(N(3),y) 0 diff(N(4),y);
diff(N(1),y) diff(N(1),x) diff(N(2),y) diff(N(2),x) diff(N(3),y) diff(N(3),x) diff(N(4),y) diff(N(4),x)];
% n5 = 1/2 * (1 - (x/lx)^2) * (1 - y/ly);
% n6 = 1/2 * (1 + (x/lx)) * (1 - (y/ly)^2);
% n7 = 1/2 * (1 - (x/lx)^2) * (1 + (y/ly));
% n8 = 1/2 * (1 - (x/lx)) * (1 - (y/ly)^2);
% n1 = N(1) - 1/2*n5 - 1/2*n8;
% n2 = N(2) - 1/2*n5 - 1/2*n6;
% n3 = N(3) - 1/2*n6 - 1/2*n7;
% n4 = N(4) - 1/2*n7 - 1/2*n8;
% % Be= [diff(n1,x), 0 , diff(n2,x), 0 , diff(n3,x), 0 , diff(n4,x), 0 , diff(n5,x), 0 , diff(n6,x), 0 , diff(n7,x), 0 , diff(n8,x), 0;
% % 0 , diff(n1,y), 0 , diff(n2,y), 0 , diff(n3,y), 0 , diff(n4,y), 0 , diff(n5,y), 0 , diff(n6,y), 0 , diff(n7,y), 0 , diff(n8,y);
% % diff(n1,y), diff(n1,x), diff(n2,y), diff(n2,x), diff(n3,y), diff(n3,x), diff(n4,y), diff(n4,x), diff(n5,y), diff(n5,x), diff(n6,y), diff(n6,x), diff(n7,y), diff(n7,x), diff(n8,y), diff(n8,x)];
% Be= [diff(n1,x), 0 , diff(n5,x), 0 , diff(n2,x), 0 , diff(n6,x), 0 , diff(n3,x), 0 , diff(n7,x), 0 , diff(n4,x), 0 , diff(n8,x), 0;
% 0 , diff(n1,y), 0 , diff(n5,y), 0 , diff(n2,y), 0 , diff(n6,y), 0 , diff(n3,y), 0 , diff(n7,y), 0 , diff(n4,y), 0 , diff(n8,y);
% diff(n1,y), diff(n1,x), diff(n5,y), diff(n5,x), diff(n2,y), diff(n2,x), diff(n6,y), diff(n6,x), diff(n3,y), diff(n3,x), diff(n7,y), diff(n7,x), diff(n4,y), diff(n4,x), diff(n8,y), diff(n8,x)];
% syms E mu
% syms d11 d12 d13 d22 d23 d33
% D = [d11 d12 d13
% d12 d22 d23
% d13 d23 d33];
E=1; mu=0.3; nu = 0.3;
D=(E/(1-mu^2))*[1 mu 0
mu 1 0
0 0 (1-mu)/2]; % E = 0.0387
% syms u1 u2 u3 u4 u5 u6 u7 u8
% u = [u1 u2 u3 u4 u5 u6 u7 u8]
% lx = 1; ly = 1;
syms x y lx ly
N(1)=(1-x/lx)*(1-y/ly);
N(2)=x/lx*(1-y/ly);
N(3)=x/lx*y/ly;
N(4)=(1-x/lx)*y/ly;
Be= [diff(N(1),x) 0 diff(N(2),x) 0 diff(N(3),x) 0 diff(N(4),x) 0;
0 diff(N(1),y) 0 diff(N(2),y) 0 diff(N(3),y) 0 diff(N(4),y);
diff(N(1),y) diff(N(1),x) diff(N(2),y) diff(N(2),x) diff(N(3),y) diff(N(3),x) diff(N(4),y) diff(N(4),x)];
% K=Be.'*D*Be;
% Ke=sym(zeros(8,8));
K=Be;
He = sym(zeros(3,8));
for i=1:3
for j=1:8
He(i,j) = int(int(K(i,j),'x',0.0,1),'y',0.0,1); % changing
% Ke(i,j)=int(int(K(i,j),'x',[0,lx]),'y',[0,ly]);
end
end
He
vpa(Ke,4)
% B divide 2
Se = ...
[- d11/2 - d13/4, - d12/2 - d13/4, d11/2 - d13/4, d13/4 - d12/2, d11/2 + d13/4, d12/2 + d13/4, d13/4 - d11/2, d12/2 - d13/4
- d12/2 - d23/4, - d22/2 - d23/4, d12/2 - d23/4, d23/4 - d22/2, d12/2 + d23/4, d22/2 + d23/4, d23/4 - d12/2, d22/2 - d23/4
- d13/2 - d33/4, - d23/2 - d33/4, d13/2 - d33/4, d33/4 - d23/2, d13/2 + d33/4, d23/2 + d33/4, d33/4 - d13/2, d23/2 - d33/4;];
% a = 1; b = 1;
% temp = E/(a*b*(1-nu^2))*[1/3*(b^2+(1-nu)/2*a^2), a*b/8*(1+nu), -1/3*(b^2-(1-nu)/4*a^2), -a*b/8*(1-3*nu), -1/6*(b^2+(1-nu)/2*a^2), -a*b/8*(1+nu), 1/6*(b^2-(1-nu)*a^2), a*b/8*(1-3*nu);
% a*b/8*(1+nu), 1/3*(a^2+(1-nu)/2*b^2), a*b/8*(1-3*nu), -1/6*(a^2-(1-nu)*b^2), -a*b/8*(1+nu), -1/6*(a^2+(1-nu)/2*b^2), -a*b/8*(1-3*nu), 1/3*(-a^2+(1-nu)/4*b^2);
% -1/3*(b^2-(1-nu)/4*a^2), a*b/8*(1-3*nu), 1/3*(b^2+(1-nu)/2*a^2), -a*b/8*(1+nu), 1/6*(b^2-(1-nu)*a^2), -a*b/8*(1-3*nu), -1/6*(b^2+(1-nu)/2*a^2), a*b/8*(1+nu);
% -a*b/8*(1-3*nu), -1/6*(a^2-(1-nu)*b^2), -a*b/8*(1+nu), 1/3*(a^2+(1-nu)/2*b^2), a*b/8*(1-3*nu), 1/3*(-a^2+(1-nu)/4*b^2), a*b/8*(1+nu), -1/6*(a^2+(1-nu)/2*b^2);
% -1/6*(b^2+(1-nu)/2*a^2), -a*b/8*(1+nu), 1/6*(b^2-(1-nu)*a^2), a*b/8*(1-3*nu), 1/3*(b^2+(1-nu)/2*a^2), a*b/8*(1+nu), 1/3*(-b^2+(1-nu)/4*a^2), -a*b/8*(1-3*nu);
% -a*b/8*(1+nu), -1/6*(a^2+(1-nu)/2*b^2), -a*b/8*(1-3*nu), 1/3*(-a^2+(1-nu)/4*b^2), a*b/8*(1+nu), 1/3*(a^2+(1-nu)/2*b^2), a*b/8*(1-3*nu), 1/6*(a^2-(1-nu)*b^2);
% 1/6*(b^2-(1-nu)*a^2), -a*b/8*(1-3*nu), -1/6*(b^2+(1-nu)/2*a^2), a*b/8*(1+nu), 1/3*(-b^2+(1-nu)/4*a^2), a*b/8*(1-3*nu), 1/3*(b^2+(1-nu)/2*a^2), -a*b/8*(1+nu);
% a*b/8*(1-3*nu), 1/3*(-a^2+(1-nu)/4*b^2), a*b/8*(1+nu), -1/6*(a^2+(1-nu)/2*b^2), -a*b/8*(1-3*nu), 1/6*(a^2-(1-nu)*b^2), -a*b/8*(1+nu), 1/3*(a^2+(1-nu)/2*b^2);]
% K = vpa(Ke,4);
% B divide 2
K = ...
[ 0.3333*d11 + 0.25*d13 + 0.08333*d33, 0.25*d12 + 0.1667*d13 + 0.1667*d23 + 0.0625*d33, 0.04167*d33 - 0.3333*d11, 0.25*d12 - 0.1667*d13 + 0.08333*d23 - 0.0625*d33, - 0.1667*d11 - 0.25*d13 - 0.04167*d33, - 0.25*d12 - 0.08333*d13 - 0.08333*d23 - 0.0625*d33, 0.1667*d11 - 0.08333*d33, 0.08333*d13 - 0.25*d12 - 0.1667*d23 + 0.0625*d33;
0.25*d12 + 0.1667*d13 + 0.1667*d23 + 0.0625*d33, 0.3333*d22 + 0.25*d23 + 0.08333*d33, 0.08333*d23 - 0.1667*d13 - 0.25*d12 + 0.0625*d33, 0.1667*d22 - 0.08333*d33, - 0.25*d12 - 0.08333*d13 - 0.08333*d23 - 0.0625*d33, - 0.1667*d22 - 0.25*d23 - 0.04167*d33, 0.25*d12 + 0.08333*d13 - 0.1667*d23 - 0.0625*d33, 0.04167*d33 - 0.3333*d22;
0.04167*d33 - 0.3333*d11, 0.08333*d23 - 0.1667*d13 - 0.25*d12 + 0.0625*d33, 0.3333*d11 - 0.25*d13 + 0.08333*d33, 0.1667*d13 - 0.25*d12 + 0.1667*d23 - 0.0625*d33, 0.1667*d11 - 0.08333*d33, 0.25*d12 + 0.08333*d13 - 0.1667*d23 - 0.0625*d33, 0.25*d13 - 0.1667*d11 - 0.04167*d33, 0.25*d12 - 0.08333*d13 - 0.08333*d23 + 0.0625*d33;
0.25*d12 - 0.1667*d13 + 0.08333*d23 - 0.0625*d33, 0.1667*d22 - 0.08333*d33, 0.1667*d13 - 0.25*d12 + 0.1667*d23 - 0.0625*d33, 0.3333*d22 - 0.25*d23 + 0.08333*d33, 0.08333*d13 - 0.25*d12 - 0.1667*d23 + 0.0625*d33, 0.04167*d33 - 0.3333*d22, 0.25*d12 - 0.08333*d13 - 0.08333*d23 + 0.0625*d33, 0.25*d23 - 0.1667*d22 - 0.04167*d33;
- 0.1667*d11 - 0.25*d13 - 0.04167*d33, - 0.25*d12 - 0.08333*d13 - 0.08333*d23 - 0.0625*d33, 0.1667*d11 - 0.08333*d33, 0.08333*d13 - 0.25*d12 - 0.1667*d23 + 0.0625*d33, 0.3333*d11 + 0.25*d13 + 0.08333*d33, 0.25*d12 + 0.1667*d13 + 0.1667*d23 + 0.0625*d33, 0.04167*d33 - 0.3333*d11, 0.25*d12 - 0.1667*d13 + 0.08333*d23 - 0.0625*d33;
-0.25*d12 - 0.08333*d13 - 0.08333*d23 - 0.0625*d33, - 0.1667*d22 - 0.25*d23 - 0.04167*d33, 0.25*d12 + 0.08333*d13 - 0.1667*d23 - 0.0625*d33, 0.04167*d33 - 0.3333*d22, 0.25*d12 + 0.1667*d13 + 0.1667*d23 + 0.0625*d33, 0.3333*d22 + 0.25*d23 + 0.08333*d33, 0.08333*d23 - 0.1667*d13 - 0.25*d12 + 0.0625*d33, 0.1667*d22 - 0.08333*d33;
0.1667*d11 - 0.08333*d33, 0.25*d12 + 0.08333*d13 - 0.1667*d23 - 0.0625*d33, 0.25*d13 - 0.1667*d11 - 0.04167*d33, 0.25*d12 - 0.08333*d13 - 0.08333*d23 + 0.0625*d33, 0.04167*d33 - 0.3333*d11, 0.08333*d23 - 0.1667*d13 - 0.25*d12 + 0.0625*d33, 0.3333*d11 - 0.25*d13 + 0.08333*d33, 0.1667*d13 - 0.25*d12 + 0.1667*d23 - 0.0625*d33;
0.08333*d13 - 0.25*d12 - 0.1667*d23 + 0.0625*d33, 0.04167*d33 - 0.3333*d22, 0.25*d12 - 0.08333*d13 - 0.08333*d23 + 0.0625*d33, 0.25*d23 - 0.1667*d22 - 0.04167*d33, 0.25*d12 - 0.1667*d13 + 0.08333*d23 - 0.0625*d33, 0.1667*d22 - 0.08333*d33, 0.1667*d13 - 0.25*d12 + 0.1667*d23 - 0.0625*d33, 0.3333*d22 - 0.25*d23 + 0.08333*d33;];
% B not divide 2
K = ...
[ d11/3 + d13/2 + d33/3, d12/4 + d13/3 + d23/3 + d33/4, d33/6 - d11/3, d12/4 - d13/3 + d23/6 - d33/4, - d11/6 - d13/2 - d33/6, - d12/4 - d13/6 - d23/6 - d33/4, d11/6 - d33/3, d13/6 - d12/4 - d23/3 + d33/4;
d12/4 + d13/3 + d23/3 + d33/4, d22/3 + d23/2 + d33/3, d23/6 - d13/3 - d12/4 + d33/4, d22/6 - d33/3, - d12/4 - d13/6 - d23/6 - d33/4, - d22/6 - d23/2 - d33/6, d12/4 + d13/6 - d23/3 - d33/4, d33/6 - d22/3;
d33/6 - d11/3, d23/6 - d13/3 - d12/4 + d33/4, d11/3 - d13/2 + d33/3, d13/3 - d12/4 + d23/3 - d33/4, d11/6 - d33/3, d12/4 + d13/6 - d23/3 - d33/4, d13/2 - d11/6 - d33/6, d12/4 - d13/6 - d23/6 + d33/4;
d12/4 - d13/3 + d23/6 - d33/4, d22/6 - d33/3, d13/3 - d12/4 + d23/3 - d33/4, d22/3 - d23/2 + d33/3, d13/6 - d12/4 - d23/3 + d33/4, d33/6 - d22/3, d12/4 - d13/6 - d23/6 + d33/4, d23/2 - d22/6 - d33/6;
- d11/6 - d13/2 - d33/6, - d12/4 - d13/6 - d23/6 - d33/4, d11/6 - d33/3, d13/6 - d12/4 - d23/3 + d33/4, d11/3 + d13/2 + d33/3, d12/4 + d13/3 + d23/3 + d33/4, d33/6 - d11/3, d12/4 - d13/3 + d23/6 - d33/4;
- d12/4 - d13/6 - d23/6 - d33/4, - d22/6 - d23/2 - d33/6, d12/4 + d13/6 - d23/3 - d33/4, d33/6 - d22/3, d12/4 + d13/3 + d23/3 + d33/4, d22/3 + d23/2 + d33/3, d23/6 - d13/3 - d12/4 + d33/4, d22/6 - d33/3;
d11/6 - d33/3, d12/4 + d13/6 - d23/3 - d33/4, d13/2 - d11/6 - d33/6, d12/4 - d13/6 - d23/6 + d33/4, d33/6 - d11/3, d23/6 - d13/3 - d12/4 + d33/4, d11/3 - d13/2 + d33/3, d13/3 - d12/4 + d23/3 - d33/4;
d13/6 - d12/4 - d23/3 + d33/4, d33/6 - d22/3, d12/4 - d13/6 - d23/6 + d33/4, d23/2 - d22/6 - d33/6, d12/4 - d13/3 + d23/6 - d33/4, d22/6 - d33/3, d13/3 - d12/4 + d23/3 - d33/4, d22/3 - d23/2 + d33/3];
% Ke = [ (E*(mu - 9))/(24*(mu^2 - 1)), -(E*(7*mu + 1))/(32*(mu^2 - 1)), (E*(mu + 15))/(48*(mu^2 - 1)), -(E*(9*mu - 1))/(32*(mu^2 - 1)), -(E*(mu - 9))/(48*(mu^2 - 1)), (E*(7*mu + 1))/(32*(mu^2 - 1)), -(E*(mu + 3))/(24*(mu^2 - 1)), (E*(9*mu - 1))/(32*(mu^2 - 1))
% -(E*(7*mu + 1))/(32*(mu^2 - 1)), (E*(mu - 9))/(24*(mu^2 - 1)), (E*(9*mu - 1))/(32*(mu^2 - 1)), -(E*(mu + 3))/(24*(mu^2 - 1)), (E*(7*mu + 1))/(32*(mu^2 - 1)), -(E*(mu - 9))/(48*(mu^2 - 1)), -(E*(9*mu - 1))/(32*(mu^2 - 1)), (E*(mu + 15))/(48*(mu^2 - 1))
% (E*(mu + 15))/(48*(mu^2 - 1)), (E*(9*mu - 1))/(32*(mu^2 - 1)), (E*(mu - 9))/(24*(mu^2 - 1)), (E*(7*mu + 1))/(32*(mu^2 - 1)), -(E*(mu + 3))/(24*(mu^2 - 1)), -(E*(9*mu - 1))/(32*(mu^2 - 1)), -(E*(mu - 9))/(48*(mu^2 - 1)), -(E*(7*mu + 1))/(32*(mu^2 - 1))
% -(E*(9*mu - 1))/(32*(mu^2 - 1)), -(E*(mu + 3))/(24*(mu^2 - 1)), (E*(7*mu + 1))/(32*(mu^2 - 1)), (E*(mu - 9))/(24*(mu^2 - 1)), (E*(9*mu - 1))/(32*(mu^2 - 1)), (E*(mu + 15))/(48*(mu^2 - 1)), -(E*(7*mu + 1))/(32*(mu^2 - 1)), -(E*(mu - 9))/(48*(mu^2 - 1))
% -(E*(mu - 9))/(48*(mu^2 - 1)), (E*(7*mu + 1))/(32*(mu^2 - 1)), -(E*(mu + 3))/(24*(mu^2 - 1)), (E*(9*mu - 1))/(32*(mu^2 - 1)), (E*(mu - 9))/(24*(mu^2 - 1)), -(E*(7*mu + 1))/(32*(mu^2 - 1)), (E*(mu + 15))/(48*(mu^2 - 1)), -(E*(9*mu - 1))/(32*(mu^2 - 1))
% (E*(7*mu + 1))/(32*(mu^2 - 1)), -(E*(mu - 9))/(48*(mu^2 - 1)), -(E*(9*mu - 1))/(32*(mu^2 - 1)), (E*(mu + 15))/(48*(mu^2 - 1)), -(E*(7*mu + 1))/(32*(mu^2 - 1)), (E*(mu - 9))/(24*(mu^2 - 1)), (E*(9*mu - 1))/(32*(mu^2 - 1)), -(E*(mu + 3))/(24*(mu^2 - 1))
% -(E*(mu + 3))/(24*(mu^2 - 1)), -(E*(9*mu - 1))/(32*(mu^2 - 1)), -(E*(mu - 9))/(48*(mu^2 - 1)), -(E*(7*mu + 1))/(32*(mu^2 - 1)), (E*(mu + 15))/(48*(mu^2 - 1)), (E*(9*mu - 1))/(32*(mu^2 - 1)), (E*(mu - 9))/(24*(mu^2 - 1)), (E*(7*mu + 1))/(32*(mu^2 - 1))
% (E*(9*mu - 1))/(32*(mu^2 - 1)), (E*(mu + 15))/(48*(mu^2 - 1)), -(E*(7*mu + 1))/(32*(mu^2 - 1)), -(E*(mu - 9))/(48*(mu^2 - 1)), -(E*(9*mu - 1))/(32*(mu^2 - 1)), -(E*(mu + 3))/(24*(mu^2 - 1)), (E*(7*mu + 1))/(32*(mu^2 - 1)), (E*(mu - 9))/(24*(mu^2 - 1))];
% Ke
% syms E mu
% ke with B is original (not divide 2)
% ke = [(E*(mu - 3)) / (6 * (mu*mu - 1)), -E / (8 * (mu - 1)), (E*(mu + 3)) / (12 * (mu*mu - 1)), -(E*(3 * mu - 1)) / (8 * (mu*mu - 1)), -(E*(mu - 3)) / (12 * (mu*mu - 1)), E / (8 * (mu - 1)), -(E*mu) / (6 * (mu*mu - 1)), (E*(3 * mu - 1)) / (8 * (mu*mu - 1))
% -E / (8 * (mu - 1)), (E*(mu - 3)) / (6 * (mu*mu - 1)), (E*(3 * mu - 1)) / (8 * (mu*mu - 1)), -(E*mu) / (6 * (mu*mu - 1)), E / (8 * (mu - 1)), -(E*(mu - 3)) / (12 * (mu*mu - 1)), -(E*(3 * mu - 1)) / (8 * (mu*mu - 1)), (E*(mu + 3)) / (12 * (mu*mu - 1))
% (E*(mu + 3)) / (12 * (mu*mu - 1)), (E*(3 * mu - 1)) / (8 * (mu*mu - 1)), (E*(mu - 3)) / (6 * (mu*mu - 1)), E / (8 * (mu - 1)), -(E*mu) / (6 * (mu*mu - 1)), -(E*(3 * mu - 1)) / (8 * (mu*mu - 1)), -(E*(mu - 3)) / (12 * (mu*mu - 1)), -E / (8 * (mu - 1))
% -(E*(3 * mu - 1)) / (8 * (mu*mu - 1)), -(E*mu) / (6 * (mu*mu - 1)), E / (8 * (mu - 1)), (E*(mu - 3)) / (6 * (mu*mu - 1)), (E*(3 * mu - 1)) / (8 * (mu*mu - 1)), (E*(mu + 3)) / (12 * (mu*mu - 1)), -E / (8 * (mu - 1)), -(E*(mu - 3)) / (12 * (mu*mu - 1))
% -(E*(mu - 3)) / (12 * (mu*mu - 1)), E / (8 * (mu - 1)), -(E*mu) / (6 * (mu*mu - 1)), (E*(3 * mu - 1)) / (8 * (mu*mu - 1)), (E*(mu - 3)) / (6 * (mu*mu - 1)), -E / (8 * (mu - 1)), (E*(mu + 3)) / (12 * (mu*mu - 1)), -(E*(3 * mu - 1)) / (8 * (mu*mu - 1))
% E / (8 * (mu - 1)), -(E*(mu - 3)) / (12 * (mu*mu - 1)), -(E*(3 * mu - 1)) / (8 * (mu*mu - 1)), (E*(mu + 3)) / (12 * (mu*mu - 1)), -E / (8 * (mu - 1)), (E*(mu - 3)) / (6 * (mu*mu - 1)), (E*(3 * mu - 1)) / (8 * (mu*mu - 1)), -(E*mu) / (6 * (mu*mu - 1))
% -(E*mu) / (6 * (mu*mu - 1)), -(E*(3 * mu - 1)) / (8 * (mu*mu - 1)), -(E*(mu - 3)) / (12 * (mu*mu - 1)), -E / (8 * (mu - 1)), (E*(mu + 3)) / (12 * (mu*mu - 1)), (E*(3 * mu - 1)) / (8 * (mu*mu - 1)), (E*(mu - 3)) / (6 * (mu*mu - 1)), E / (8 * (mu - 1))
% (E*(3 * mu - 1)) / (8 * (mu*mu - 1)), (E*(mu + 3)) / (12 * (mu*mu - 1)), -E / (8 * (mu - 1)), -(E*(mu - 3)) / (12 * (mu*mu - 1)), -(E*(3 * mu - 1)) / (8 * (mu*mu - 1)), -(E*mu) / (6 * (mu*mu - 1)), E / (8 * (mu - 1)), (E*(mu - 3)) / (6 * (mu*mu - 1)) ]