115 lines
3.6 KiB
Matlab
115 lines
3.6 KiB
Matlab
function [r] = GridProp_pre(dx, dy,dz,nx,ny,nz,NTG)
|
|
% coordinates xy cordinate，即每个节点的x,y坐标
|
|
% nodes connection，即每个子区域对应的八个节点数
|
|
% nP numbers of point，即节点数量
|
|
% nE numbers of element，即子区域数量
|
|
nP = (nx + 1) * (ny + 1)*(nz+1);%计算节点数量
|
|
nE = nx * ny*nz;%子区域数量计算
|
|
% dx=dx*ones(nx,1);%此处是为了暂时妥协
|
|
% dy=dy*ones(nx,1);
|
|
% dz=dz*ones(nx,1);
|
|
dxv = zeros(nE, 1);%每个子区域x方向上尺寸矩阵
|
|
dyv = zeros(nE, 1);%每个子区域y方向上尺寸矩阵
|
|
dzv = zeros(nE, 1);
|
|
zm=zeros(nE, 1);%基质网格计算深度（以下平面为基准计算得到的高度，值为正数）
|
|
xm=zeros(nE, 1);
|
|
ym=zeros(nE, 1);
|
|
%给每个子区域的x方向和y方向赋尺寸
|
|
vm=zeros(nE, 1);%基质网格体积
|
|
for k=1:nz
|
|
for i = 1 : nx
|
|
for j = 1 : ny
|
|
dxv(i+(j-1)*nx+nx*ny*(k-1)) = dx(i);
|
|
dyv(i+(j-1)*nx+nx*ny*(k-1)) = dy(j);
|
|
dzv(i+(j-1)*nx+nx*ny*(k-1))= dz(k);
|
|
vm(i+(j-1)*nx+nx*ny*(k-1))=dx(i)*dy(j)*dz(k)*NTG(i+(j-1)*nx+nx*ny*(k-1));
|
|
end
|
|
end
|
|
end
|
|
coordinates = zeros(nP, 3);%初始化每个节点的x,y坐标矩阵
|
|
nodes = zeros(nE, 8);%初始化每个子区域对应的八个节点数的矩阵
|
|
%计算每个节点的坐标
|
|
x0 = zeros(nx+1,1);
|
|
y0 = zeros(ny+1,1);
|
|
z0=zeros(nz+1,1);
|
|
sumx = 0;
|
|
sumy = 0;
|
|
sumz=0;
|
|
for i = 1 : nx
|
|
sumx = sumx+dx(i);
|
|
x0(i+1) = sumx;
|
|
end
|
|
for i = 1 : ny
|
|
sumy = sumy+dy(i);
|
|
y0(i+1) = sumy;
|
|
end
|
|
for i = 1 : nz
|
|
sumz = sumz+dz(i);
|
|
z0(i+1) = sumz;
|
|
end
|
|
for k=1:nz+1
|
|
for j = 1 : ny + 1
|
|
for i = 1 : nx + 1
|
|
coordinates(i + (j - 1) * (nx + 1)+(nx+1)*(ny+1)*(k-1), 1) = x0(i);%(i,j）节点的x坐标赋值
|
|
coordinates(i + (j - 1) * (nx + 1)+(nx+1)*(ny+1)*(k-1), 2) = y0(j);%(i,j）节点的y坐标赋值
|
|
coordinates(i + (j - 1) * (nx + 1)+(nx+1)*(ny+1)*(k-1), 3) = z0(k);
|
|
end
|
|
end
|
|
end
|
|
%计算每个子区域的八个顶点(节点）对应的节点编号
|
|
for k=1:nz
|
|
for j = 1 : ny
|
|
for i = 1 : nx
|
|
nodes(i + (j - 1) * nx+nx*ny*(k-1), 1) = i + (j - 1) * (nx + 1)+(k-1)*(nx+1)*(ny+1);
|
|
nodes(i + (j - 1) * nx+nx*ny*(k-1), 2) = nodes(i + (j - 1) * nx+nx*ny*(k-1), 1) + 1;
|
|
nodes(i + (j - 1) * nx+nx*ny*(k-1), 3) = nodes(i + (j - 1) * nx+nx*ny*(k-1), 1) + nx + 1;
|
|
nodes(i + (j - 1) * nx+nx*ny*(k-1), 4) = nodes(i + (j - 1) * nx+nx*ny*(k-1), 1) + nx + 2;
|
|
nodes(i + (j - 1) * nx+nx*ny*(k-1), 5) = nodes(i + (j - 1) * nx+nx*ny*(k-1), 1)+(nx+1)*(ny+1);
|
|
nodes(i + (j - 1) * nx+nx*ny*(k-1), 6) = nodes(i + (j - 1) * nx+nx*ny*(k-1), 5) + 1;
|
|
nodes(i + (j - 1) * nx+nx*ny*(k-1), 7) = nodes(i + (j - 1) * nx+nx*ny*(k-1), 5) + nx + 1;
|
|
nodes(i + (j - 1) * nx+nx*ny*(k-1), 8) = nodes(i + (j - 1) * nx+nx*ny*(k-1), 5) + nx + 2;
|
|
zm(i+(j-1)*nx+nx*ny*(k-1))=sum(coordinates(nodes(i + (j - 1) * nx+nx*ny*(k-1), :),3))/8;
|
|
xm(i+(j-1)*nx+nx*ny*(k-1))=sum(coordinates(nodes(i + (j - 1) * nx+nx*ny*(k-1), :),1))/8;
|
|
ym(i+(j-1)*nx+nx*ny*(k-1))=sum(coordinates(nodes(i + (j - 1) * nx+nx*ny*(k-1), :),2))/8;
|
|
end
|
|
end
|
|
end
|
|
|
|
cell_mid_coordinates = [xm,ym,zm];
|
|
|
|
xrao=[0];
|
|
for i=1:nx
|
|
xxiang=sum(dx(1,1:i));
|
|
xrao=[xrao xxiang];
|
|
end
|
|
yrao=[0];
|
|
for i=1:ny
|
|
yxiang=sum(dy(1,1:i));
|
|
yrao=[yrao yxiang];
|
|
end
|
|
zrao=[0];
|
|
for i=1:nz
|
|
zxiang=sum(dz(1,1:i));
|
|
zrao=[zrao zxiang];
|
|
end
|
|
r.dx = dx;
|
|
r.dy = dy;
|
|
r.dz = dz;
|
|
r.nx = nx;
|
|
r.ny = ny;
|
|
r.nz = nz;
|
|
r.coordinates = coordinates;
|
|
r.nodes = nodes;
|
|
r.nP = nP;
|
|
r.nmc = nE;
|
|
r.dxv = dxv;
|
|
r.dyv = dyv;
|
|
r.dzv = dzv;
|
|
r.zm = zm;
|
|
r.vm = vm;
|
|
r.xrao = xrao;
|
|
r.yrao = yrao;
|
|
r.zrao = zrao;
|
|
r.cell_mid_coords = cell_mid_coordinates;
|
|
r.NTG = NTG;
|
|
end |