Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Welcome To Ask or Share your Answers For Others

Categories

0 votes
782 views
in Technique[技术] by (71.8m points)

calculate 3d strain using systems of linear equation (MATLAB)

I have 10 elements with their serial number and x,y,z coordinates in their initial and final position : data0 and data1 respectively.

I intend to calculate the strain as follows:

  1. for each element in data0, I will search for 3 elements closest to it and for each of the four closest elements (i, j, k, l), the displacements ui, uj, uk, ul in the x direction; vi, vj, vk, vl in the y direction; and wi, wj, wk, wl in the z direction are all calculated

  2. simultaneous equations are created for each 4 element set to obtain a0, a1, a2; b0 ,b1 ,b2; c0 ,c1 ,c2
    "-[1, xi,yi,zi;1,xj,yj,zj;1,xk,yk,zk;1,xl,yl,zl] * [a0,a1,a2,a3] = [ui,uj,uk,ul] -[1, xi,yi,zi;1,xj,yj,zj;1,xk,yk,zk;1,xl,yl,zl] * [b0,b1,b2,b3] = [vi,vj,vk,vl] -[1, xi,yi,zi;1,xj,yj,zj;1,xk,yk,zk;1,xl,yl,zl] * [c0,c1,c2,c3] = [wi,wj,wk,wl]"

  3. finally for each element the following strain is calculated using the equation in the calculate micro strain section of the code below.

I tried creating the following code but got my simultaneous equation loop is not well structured

data0 = [1, 46.15, 29.55, 4.36; 
2, 56.19, 34.22, 5.33;
3, 39.01, 33.28, 7.27; 
4, 29.51, 37.41, 5.89;
5, 65.51, 37.97, 5.04;
6, 50.78, 47.22, 3.25; 
7, 72.69, 47.58, 4.34; 
8, 30.07, 48.18, 5.62; 
9, 41.34, 50.11, 3.91;
10,61.13, 52.16, 2.22]                                                                                                                                                                                                                                  
%Extract element number and store element_number
Element_No1=data0(:,1);
X1=data0(:,2);
Y1=data0(:,3);
Z1=data0(:,4);
EUC_DIST= data0(:,5);

data1= [1,46.83,30,4.69;
2,56.84,35.49,6.51;
3,39.39,33.56,7.11;
4,29.74,38.16,5.87;
5,64.5,38.55,5.61;
6,50.74,47.57,3.22;
7,71.88,47.8,4.47;
8,29.72,48.2,5.59;
9,41.13,50.23,3.88;
10,60.61,52.4,2.29]

Element_No2=data1(:,1);
X2=data1(:,2);
Y2=data1(:,3);
Z2=data1(:,4);
%Dx= data1(:,5);
%Dy= data1(:,6);
%Dz= data1(:,7);

%Check matrix size           
[m,n]=size(data0);
[p,o]=size(data1);

%count the number of data points
Q0 = numel(data0(:,1));
Q1 = numel(data1(:,1));
q1 = [1:Q1]';

%create matrix of distance difference with zeros
DIST_DIFF = zeros(Q0,Q1);

%Calculate the eucledean distsnce difference between each element and others
for j=1:m
for i=1:m
DIST_DIFF(j,i) = abs((data0(j,5)-data0(i,5)));
DIST_DIFF(DIST_DIFF == 0) = NaN;
end 
end

%add the element number to the calculated distance difference 
DIST_DIFF = horzcat(DIST_DIFF,q1);


%Get the closest 3 element counterpart and store them in EUCL_DIST
[r,s]=size(DIST_DIFF);
A = [DIST_DIFF(:,1),DIST_DIFF(:,(Q0+1))];
SA = size(A);
NA = zeros(SA);
EUCL_DIST = [];

for i= 1:s-1
%generate matrix for each column of DIST_DIFF and Attach element number to each generated column
NA = [DIST_DIFF(:,i),DIST_DIFF(:,(Q0+1))];
%sort each generated column 
S_NA= sortrows(NA);
%take out the 3 closest element
T= [S_NA(1:3,1:2)]; %3 here indicate number of elements to display 
EUCL_DIST = [EUCL_DIST,T];
end


%Attach other element properties (x,y,z) to matched particles ID

%Isolate the 3 closest particle number and leave out the calculated distance
y0= EUCL_DIST(:,2:2:end);

[r,s]=size(y0);
TOG = [];

for i= 1:s
%generate matrix for each column of DIST_DIFF and Attach element number to each generated column 
%take out the 3 closest element
q=ismember(data0(:,1), y0(:,i), 'rows');
z=find(q);
Close_Elements =data0(z,1:4);
TOG = [TOG,Close_Elements];
end


ALL_VAR = zeros(Q0,21);

[p,t]=size(TOG);
[g,h]=size(data0);
[q,r]=size(data1);

for i= 1:t
for ii= 1:h
for iii= 1:r
A1 =[1, data0(1,2), data0(1,3), data0(1,4);
1, TOG(1,2), TOG(1,3),  TOG(1,4);
1, TOG(2,2), TOG(2,3),  TOG(2,4);
1, TOG(3,2), TOG(3,3),  TOG(3,4)]; 
        
A2 =[1, data0(1,2), data0(1,3), data0(1,4);
1,  TOG(1,2), TOG(1,3),TOG(1,4);
1, TOG(2,2),TOG(2,3),TOG(2,4);
1, TOG(3,2), TOG(3,3),TOG(3,4)]; 

A3 =[1, data0(1,2), data0(1,3), data0(1,4);
1,  TOG(1,2),TOG(1,3),TOG(1,4);
1, TOG(2,2),TOG(2,3),TOG(2,4);
1, TOG(3,2),TOG(3,3),TOG(3,4)]; 


B1 = [(data0(1,2)-data1(1,2));
(data0(2,2)-data1(2,2));
(data0(3,2)-data1(3,2));   
(data0(4,2)-data1(4,2))];
        
B2 = [(data0(1,3)-data1(1,3));
(data0(2,3)-data1(2,3));
(data0(3,3)-data1(3,3));
(data0(4,3)-data1(4,3))];
        
B3 = [(data0(1,4)-data1(1,4));
(data0(2,4)-data1(2,4));
(data0(3,4)-data1(3,4));
(data0(4,4)-data1(4,4))];

a = linsolve(A1,B1);
b = linsolve(A2,B2);
c = linsolve(A3,B3);

Const = [a,b,c];
end
end
end
%%%%%Difference between initial and final particle position
Du = X2 - X1;
Dv= y2 - y1;
Dw= z2 - z1;
%%%Difference between each close particle and the particle under consideration
Dx= x2nd -x1;
Dy= x3rd-x1;
Dz= x4th-x1;

%Calculate micro strain - update du ,dv ,dw ,dx ,dy ,dz to ALL_VARIABLES
a1 = du/dx ;
b2 = dv/dy;
c3 = dw/dz;

Ex = a1;
Ey = b2;
Ez = c3;

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome To Ask or Share your Answers For Others

1 Reply

0 votes
by (71.8m points)
等待大神答复

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
OGeek|极客中国-欢迎来到极客的世界,一个免费开放的程序员编程交流平台!开放,进步,分享!让技术改变生活,让极客改变未来! Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Click Here to Ask a Question

...