-
Notifications
You must be signed in to change notification settings - Fork 56
/
cellBoundary3D.m
executable file
·123 lines (112 loc) · 2.98 KB
/
cellBoundary3D.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
function phiBC = cellBoundary3D(phi, BC)
% function phiBC = cellBoundary2D(phi, BC)
% It creates the matrix of coefficient based on the BC structure provided
% by the user. It also generates the right hand side vector of the linear
% system of equations
%
% SYNOPSIS:
% phiBC = cellBoundary2D(phi, BC)
%
% PARAMETERS:
% BC: boundary condition structure created by createBC function
% phi: cell variable created by createCellVariable
%
% RETURNS:
% phiBC: a cell variable including the values of the ghost cells
%
% EXAMPLE:
% m = createMesh2D(3,4,1,1);
% phi = createCellVariable(m,1);
% bc = createBC(m);
% phi_with_ghost = cellBoundary(phi, bc)
%
% SEE ALSO:
%
% extract data from the mesh structure
Nxyz = BC.domain.dims;
Nx = Nxyz(1); Ny = Nxyz(2); Nz = Nxyz(3);
dx_1 = BC.domain.cellsize.x(1);
dx_end = BC.domain.cellsize.x(end);
dy_1 = BC.domain.cellsize.y(1);
dy_end = BC.domain.cellsize.y(end);
dz_1 = BC.domain.cellsize.z(1);
dz_end = BC.domain.cellsize.z(end);
% define the output matrix
phiBC = zeros(Nx+2, Ny+2, Nz+2);
phiBC(2:Nx+1, 2:Ny+1, 2:Nz+1) = phi;
% Assign values to the boundary values
if (BC.top.periodic==0) && (BC.bottom.periodic==0)
% top boundary
j=Ny+2;
i = 2:Nx+1;
k = 2:Nz+1;
phiBC(i,j,k)= ...
(BC.top.c-squeeze(phi(:,end,:)).*(-BC.top.a/dy_end+BC.top.b/2))./(BC.top.a/dy_end+BC.top.b/2);
% Bottom boundary
j=1;
i = 2:Nx+1;
k = 2:Nz+1;
phiBC(i,j,k)= ...
(BC.bottom.c-squeeze(phi(:,1,:)).*(BC.bottom.a/dy_1+BC.bottom.b/2))./(-BC.bottom.a/dy_1+BC.bottom.b/2);
else
% top boundary
j=Ny+2;
i = 2:Nx+1;
k = 2:Nz+1;
phiBC(i,j,k)= phi(:,1,:);
% Bottom boundary
j=1;
i = 2:Nx+1;
k = 2:Nz+1;
phiBC(i,j,k)= phi(:,end,:);
end
if (BC.left.periodic==0) && (BC.right.periodic==0)
% Right boundary
i = Nx+2;
j = 2:Ny+1;
k = 2:Nz+1;
phiBC(i,j,k)= ...
(BC.right.c-squeeze(phi(end,:,:)).*(-BC.right.a/dx_end+BC.right.b/2))./(BC.right.a/dx_end+BC.right.b/2);
% Left boundary
i = 1;
j = 2:Ny+1;
k = 2:Nz+1;
phiBC(i,j,k)= ...
(BC.left.c-squeeze(phi(1,:,:)).*(BC.left.a/dx_1+BC.left.b/2))./(-BC.left.a/dx_1+BC.left.b/2);
else
% Right boundary
i = Nx+2;
j = 2:Ny+1;
k = 2:Nz+1;
phiBC(i,j,k)= phi(1,:,:);
% Left boundary
i = 1;
j = 2:Ny+1;
k = 2:Nz+1;
phiBC(i,j,k)= phi(end,:,:);
end
if (BC.bottom.periodic==0) && (BC.top.periodic==0)
% front boundary
i = 2:Nx+1;
j = 2:Ny+1;
k = Nz+2;
phiBC(i,j,k)= ...
(BC.front.c-squeeze(phi(:,:,end)).*(-BC.front.a/dz_end+BC.front.b/2))./(BC.front.a/dz_end+BC.front.b/2);
% back boundary
i = 2:Nx+1;
j = 2:Ny+1;
k = 1;
phiBC(i,j,k)= ...
(BC.back.c-squeeze(phi(:,:,1)).*(BC.back.a/dz_1+BC.back.b/2))./(-BC.back.a/dz_1+BC.back.b/2);
else
% front boundary
i = 2:Nx+1;
j = 2:Ny+1;
k = Nz+2;
phiBC(i,j,k)= phi(:,:,1);
% back boundary
i = 2:Nx+1;
j = 2:Ny+1;
k = 1;
phiBC(i,j,k)= phi(:,:,end);
end