-
Notifications
You must be signed in to change notification settings - Fork 1
/
vt_circularTubeGeneration.m
227 lines (185 loc) · 10.1 KB
/
vt_circularTubeGeneration.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
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
function [simGridParam, PV_N, tubeStartArea, tubeStart, tubeEnd, totalTubeLengthInCells, currTubeSectionDiameterCells_SegmentCounter, listenerInfo, sourceInfo] = vt_circularTubeGeneration(simulationType, junctionType, vowel, boundaryInterpolation, gridCellTypes, baffleSwitchFlag, pmlSwitch, pmlLayer, rad, ds)
% Define units
meter = 1;
centimeter = 1e-2*meter;
milimeter = 1e-3*meter;
diameter_mul = 1;
% STEP XX: Define tube geometry
[segment_length, tubeLength, tubeSectionArea_incm2] = vt_areaFunction(simulationType, junctionType, vowel);
% Tube section area in m^2
tubeSectionArea_inm2 = tubeSectionArea_incm2.*(centimeter*centimeter);
% Virtual mic/listener position near the mouth-end [inside the tube]
micXdistanceFromMouthEnd = 3*milimeter;
micxdistanceInCellsFromMouthEnd = round(micXdistanceFromMouthEnd/ds);
% number of tube segments
numSections = length(tubeSectionArea_inm2);
% Calculate tube section diameters
tubeSectionOriginalDiameter = 2*sqrt(tubeSectionArea_inm2./pi).* diameter_mul;
tubeSectionDiameterCells = round(tubeSectionOriginalDiameter./ds);
% Cross-sectional area at the tube start end
% I didn't use "tubeSectionArea_inm2" directly because I wanted to use
% the derived area by using the derived diameter.
tubeStartArea = pi*(tubeSectionOriginalDiameter(1)/2)^2;
% Change the tubeSectionDiameterCells to 1 if it contains 0
tubeSectionDiameterCells(tubeSectionDiameterCells==0)=1;
%STEP XX: Choose the best possible odd number from the Diameter array
for diameterCounter = 1:numSections
% Verify if the cellsPerDiameter is odd or not
if mod(tubeSectionDiameterCells(diameterCounter), 2) == 0
% Find the difference between the rounded and the actual
% diameter value = Estimated diameter-Actual diameter
diff = tubeSectionDiameterCells(diameterCounter) - ...
(tubeSectionOriginalDiameter(diameterCounter)/ds);
if diff>0
tubeSectionDiameterCells(diameterCounter) = ...
tubeSectionDiameterCells(diameterCounter)-1;
else
tubeSectionDiameterCells(diameterCounter) = ...
tubeSectionDiameterCells(diameterCounter)+1;
end
end
end
% Number of cells for total tube length
actualTubeLength = numSections*segment_length;
totalTubeLengthInCells = round(actualTubeLength/ds);
if baffleSwitchFlag == 1
% TO DO: Deb Implement the baffle condition
% Deb: For eccentric tube condition how should we include baffle?
else
domainX = totalTubeLengthInCells + 1 + 1; % +1 is excitation and and dirichlet layer
domainY = max(tubeSectionDiameterCells) + 2; % +2 is tube walls
domainZ = max(tubeSectionDiameterCells) + 2; % +2 is tube walls
end
% Derive frame size
frameX = domainX + 2; %2 = dead cells
frameX = frameX + 2*pmlSwitch*pmlLayer;
frameY = domainY + 2; %2 = dead cells
frameY = frameY + 2*pmlSwitch*pmlLayer;
frameZ = domainZ + 2; %2 = dead cells
frameZ = frameZ + 2*pmlSwitch*pmlLayer;
% STEP XX: Add all the frame params to the simGridParam
simGridParam.frameX = frameX;
simGridParam.frameY = frameY;
simGridParam.frameZ = frameZ;
% Define PV_N array to store pressure and velocity components
% PV_N(:,:,1) = To store cell pressure
% PV_N(:,:,2) = To store Vx
% PV_N(:,:,3) = To store Vy
% PV_N(:,:,4) = To store Vz
% PV_N(:,:,5) = To store grid cell types
PV_N = zeros(frameY, frameX, frameZ, 5);
% Define cell types and store it in PV_N(,,5)
% Declare all the cells as air by default
PV_N(1:frameY, 1:frameX, 1:frameZ, 5) = gridCellTypes.cell_air;
% STEP XX: Create the regular tube geometry
% The axis of the cylindrical tube is along the x-axis, direction along
% which acoustic wave propagates. We will start building the tube model
% from the starting position towards the end position.
% STEP XX.XX Find the center of the tube along the yz-plane
centerY = ceil(frameY/2);
centerZ = ceil(frameZ/2);
% STEP XX.XX Find the starting and ending point of the tube.
% 1 is: for starting point, dead layer and excitation
% 2 is: for cell_head in case of the circular baffle and to have an
% empty cell between the cell_head and pml layers.
if rad ~=3
tubeStart.startX = 1 + 1 + pmlLayer*pmlSwitch + 1;
tubeEnd.endX = tubeStart.startX + totalTubeLengthInCells - 1 ;
else
%TO DO: Deb check this later tubeStartX = 1 + 1 + pmlLayer*pmlSwitch + 1+2;
end
tubeStart.startY = centerY;
tubeEnd.endY = centerY;
tubeStart.startZ = centerZ;
tubeEnd.endZ = centerZ;
% STEP XX: Store the cummulative length of each tube section
tubeCummSectionLength = zeros(1, numSections);
for sectionCount = 1:numSections
tubeCummSectionLength(sectionCount) = segment_length*sectionCount;
end
if junctionType == 1
[PV_N, currTubeSectionDiameterCells_SegmentCounter]...
= vt_centricCircularContourGeneration(PV_N, numSections, totalTubeLengthInCells, tubeSectionDiameterCells, tubeCummSectionLength, boundaryInterpolation, tubeStart, gridCellTypes, ds);
else
[PV_N, currTubeSectionDiameterCells_SegmentCounter]...
= vt_eccentricCircularContourGeneration(PV_N, numSections, totalTubeLengthInCells, tubeSectionDiameterCells, tubeCummSectionLength, boundaryInterpolation, tubeStart, pmlSwitch, pmlLayer, gridCellTypes, ds);
end
% STEPXX: Calculate tube midline length for eccentric geometry case
% This can be done by checking the midpoint along the y-axis
vtMidlineLength = 0;
for vtCrossSectionCounter = tubeStart.startX:tubeEnd.endX
[gridPlaneProp, gridCellTypeInplane] = vt_findCellTypes(PV_N, gridCellTypes, vtCrossSectionCounter);
airCells = find(gridPlaneProp(:, tubeStart.startZ)==gridCellTypeInplane.inVTContour);
currMidYPosition = (min(airCells) + max(airCells))/2;
if vtCrossSectionCounter == tubeStart.startX
prevMidYPosition = currMidYPosition;
continue;
end
if prevMidYPosition == currMidYPosition
vtMidlineLength = vtMidlineLength+ds;
prevMidYPosition = currMidYPosition;
else
yPosDifferenceLen = abs(currMidYPosition-prevMidYPosition)*ds;
vtMidlineLength = vtMidlineLength + sqrt(ds^2 + yPosDifferenceLen^2);
prevMidYPosition = currMidYPosition;
end
end
finalvtMidlineLength = vtMidlineLength+ds;
% Percentage error in total tube length
totalTubeLengthError = (finalvtMidlineLength-tubeLength)/...
(tubeLength);
fprintf('Approximated tube length percentage error = %.4f \n',totalTubeLengthError*100);
% STEPXX: Define excitation cells
% Check the grid cell types for the yz-plane that exists besides
% (left-side) the tube starting point
xExcitationPos = tubeStart.startX-1;
[gridPlaneProp, gridCellTypeInplane] = vt_findCellTypes(PV_N, gridCellTypes, tubeStart.startX);
% Find the grid size and traverse through yz-plane to assign
% cell_excitation
gridSize = size(PV_N);
for yCount = 1:gridSize(1)
for zCount = 1:gridSize(3)
if gridPlaneProp(yCount, zCount) == gridCellTypeInplane.inVTContour
PV_N(yCount, xExcitationPos, zCount, 5) = gridCellTypes.cell_excitation;
elseif gridPlaneProp(yCount, zCount) == gridCellTypeInplane.onVTContour
PV_N(yCount, xExcitationPos, zCount, 5) = gridCellTypes.cell_wall;
end
end
end
% STEPXX: Define no_Pressure cells [Dirichlet Boundary Condition]
if rad==2
tubeEnd.endX = tubeStart.startX + totalTubeLengthInCells - 1 ;
[gridPlaneProp, gridCellTypeInplane] = vt_findCellTypes(PV_N, gridCellTypes, tubeEnd.endX);
% Traverse through yz-plane to assign cell_noPressure
for yCount = 1:gridSize(1)
for zCount = 1:gridSize(3)
if gridPlaneProp(yCount, zCount) == gridCellTypeInplane.inVTContour
PV_N(yCount, tubeEnd.endX+1, zCount, 5) = gridCellTypes.cell_noPressure;
elseif gridPlaneProp(yCount, zCount) == gridCellTypeInplane.onVTContour
PV_N(yCount, tubeEnd.endX+1, zCount, 5) = gridCellTypes.cell_noPressure;
end
end
end
end
% Place two microphines:
% 1. Near mouth-end [3mm inside of mouth] = listenerInfo
% 2. Near glottal-end [Next to excitation cells] = sourceInfo
% Deb: We need to careful for centric and eccentric tube.
% First find the VT cross-section along the yz plane.
% Then find the y coordinate for source and listener by checking
% the min/max positions.
% STEPXX: Determine the microphone/listener position near the mouth-end
listenerInfo.listenerX = tubeEnd.endX-micxdistanceInCellsFromMouthEnd;
listenerInfo.listenerZ = tubeEnd.endZ;
[gridPlaneProp, gridCellTypeInplane] = vt_findCellTypes(PV_N, gridCellTypes, listenerInfo.listenerX);
airCells = find(gridPlaneProp(:, tubeEnd.endZ)==gridCellTypeInplane.inVTContour);
yPosition = (min(airCells) + max(airCells))/2;
listenerInfo.listenerY = yPosition;
% STEPXX: Determine the source position near the glottal-end
sourceInfo.sourceX = tubeStart.startX+1;
sourceInfo.sourceZ = tubeStart.startZ;
[gridPlaneProp, gridCellTypeInplane] = vt_findCellTypes(PV_N, gridCellTypes, sourceInfo.sourceX);
airCells = find(gridPlaneProp(:, tubeStart.startZ)==gridCellTypeInplane.inVTContour);
yPosition = (min(airCells) + max(airCells))/2;
sourceInfo.sourceY = yPosition;
end