skewWindSim/getRegularGrid.m at main · ECheynet/skewWindSim · GitHub

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
function [geometry,minZ] = getRegularGrid(minZ,Zlim,Ylim,dy,dz)

z = [];
y = [];

minZ = max(1,minZ);

if minZ-dz==0,
    warning('lowest height identified as 0 m; minZ is set at minZ = min>+1 m')
    minZ = minZ+1;
end


minY = -Ylim/2;
Zlim = Zlim+ dz;
z1 = [(minZ-dz)/2:dz:(Zlim)];
y1 = (minY-dy).*ones(1,numel(z1));

for ii=1:numel(minY:dy:-minY)
    y = [y,y1+dy*ii];
    if mod(ii,2)==1
        z = [z,fliplr(z1)];
    else
        z = [z,(z1)];
    end
end

y = y(:);
z = z(:);


x = zeros(size(y));

% Nodes
geometry.node.X = x; % must be [1 x M]
geometry.node.Y = y; % must be [1 x M]
geometry.node.Z = z;
% elements whose coordinates are at midpoint between the elements
geometry.element.X =  (geometry.node.X(2:end)+geometry.node.X(1:end-1))/2;
geometry.element.Y = (geometry.node.Y(2:end)+geometry.node.Y(1:end-1))/2;
geometry.element.Z = (geometry.node.Z(2:end)+geometry.node.Z(1:end-1))/2;


geometry.node.Z(geometry.node.Z<0) = 0.1;
geometry.element.Z(geometry.element.Z<0) = 0.1;


end