Mercurial > repos > public > sbplib

comparison +grid/Cartesian.m @ 832:5573913a0949 feature/burgers1d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merged with default, and updated +scheme/Burgers1D accordingly
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Tue, 11 Sep 2018 15:58:35 +0200
parents 031d6db97270
children
comparison
equal deleted inserted replaced
831:d0934d1143b7 832:5573913a0949
1 classdef Cartesian < grid.Structured
2 properties
3 n % Number of points in the grid
4 d % Number of dimensions
5 m % Number of points in each direction
6 x % Cell array of vectors with node placement for each dimension.
7 h % Spacing/Scaling
8 lim % Cell array of left and right boundaries for each dimension.
9 end
10
11 % General d dimensional grid with n points
12 methods
13 % Creates a cartesian grid given vectors conatining the coordinates
14 % in each direction
15 function obj = Cartesian(varargin)
16 obj.d = length(varargin);
17
18 for i = 1:obj.d
19 assert(isnumeric(varargin{i}), 'Coordinate inputs must be vectors.')
20
21 obj.x{i} = varargin{i};
22 obj.m(i) = length(varargin{i});
23 end
24
25 obj.n = prod(obj.m);
26 if obj.n == 0
27 error('grid:Cartesian:EmptyGrid','Input parameter gives an empty grid.')
28 end
29
30 obj.h = [];
31
32 obj.lim = cell(1,obj.d);
33 for i = 1:obj.d
34 obj.lim{i} = {obj.x{i}(1), obj.x{i}(end)};
35 end
36 end
37 % n returns the number of points in the grid
38 function o = N(obj)
39 o = obj.n;
40 end
41
42 % d returns the spatial dimension of the grid
43 function o = D(obj)
44 o = obj.d;
45 end
46
47 function m = size(obj)
48 m = obj.m;
49 end
50
51 % points returns a n x d matrix containing the coordianets for all points.
52 % points are ordered according to the kronecker product with X*Y*Z
53 function X = points(obj)
54 X = zeros(obj.n, obj.d);
55
56 for i = 1:obj.d
57 if iscolumn(obj.x{i})
58 c = obj.x{i};
59 else
60 c = obj.x{i}';
61 end
62
63 m_before = prod(obj.m(1:i-1));
64 m_after = prod(obj.m(i+1:end));
65
66 X(:,i) = kr(ones(m_before,1),c,ones(m_after,1));
67 end
68 end
69
70 % matrices returns a cell array with coordinates in matrix form.
71 % For 2d case these will have to be transposed to work with plotting routines.
72 function X = matrices(obj)
73
74 if obj.d == 1 % There is no 1d matrix data type in matlab, handle special case
75 X{1} = reshape(obj.x{1}, [obj.m 1]);
76 return
77 end
78
79 X = cell(1,obj.d);
80 for i = 1:obj.d
81 s = ones(1,obj.d);
82 s(i) = obj.m(i);
83
84 t = reshape(obj.x{i},s);
85
86 s = obj.m;
87 s(i) = 1;
88 X{i} = repmat(t,s);
89 end
90 end
91
92 function h = scaling(obj)
93 if isempty(obj.h)
94 error('grid:Cartesian:NoScalingSet', 'No scaling set')
95 end
96
97 h = obj.h;
98 end
99
100 % Restricts the grid function gf on obj to the subgrid g.
101 % Only works for even multiples
102 function gf = restrictFunc(obj, gf, g)
103 m1 = obj.m;
104 m2 = g.m;
105
106 % Check the input
107 if prod(m1) ~= numel(gf)
108 error('grid:Cartesian:restrictFunc:NonMatchingFunctionSize', 'The grid function has to few or too many points.');
109 end
110
111 if ~all(mod(m1-1,m2-1) == 0)
112 error('grid:Cartesian:restrictFunc:NonMatchingGrids', 'Only integer downsamplings are allowed');
113 end
114
115 % Calculate stride for each dimension
116 stride = (m1-1)./(m2-1);
117
118 % Create downsampling indecies
119 I = {};
120 for i = 1:length(m1)
121 I{i} = 1:stride(i):m1(i);
122 end
123
124 gf = reshapeRowMaj(gf, m1);
125 gf = gf(I{:});
126 gf = reshapeRowMaj(gf, prod(m2));
127 end
128
129 % Projects the grid function gf on obj to the grid g.
130 function gf = projectFunc(obj, gf, g)
131 error('grid:Cartesian:NotImplemented')
132 end
133
134 % Return the names of all boundaries in this grid.
135 function bs = getBoundaryNames(obj)
136 switch obj.D()
137 case 1
138 bs = {'l', 'r'};
139 case 2
140 bs = {'w', 'e', 's', 'n'};
141 case 3
142 bs = {'w', 'e', 's', 'n', 'd', 'u'};
143 otherwise
144 error('not implemented');
145 end
146 end
147
148 % Return coordinates for the given boundary
149 function X = getBoundary(obj, name)
150 % In what dimension is the boundary?
151 switch name
152 case {'l', 'r', 'w', 'e'}
153 D = 1;
154 case {'s', 'n'}
155 D = 2;
156 case {'d', 'u'}
157 D = 3;
158 otherwise
159 error('not implemented');
160 end
161
162 % At what index is the boundary?
163 switch name
164 case {'l', 'w', 's', 'd'}
165 index = 1;
166 case {'r', 'e', 'n', 'u'}
167 index = obj.m(D);
168 otherwise
169 error('not implemented');
170 end
171
172
173
174 I = cell(1, obj.d);
175 for i = 1:obj.d
176 if i == D
177 I{i} = index;
178 else
179 I{i} = ':';
180 end
181 end
182
183 % Calculate size of result:
184 m = obj.m;
185 m(D) = [];
186 N = prod(m);
187
188 X = zeros(N, obj.d);
189
190 coordMat = obj.matrices();
191 for i = 1:length(coordMat)
192 Xtemp = coordMat{i}(I{:});
193 X(:,i) = reshapeRowMaj(Xtemp, [N,1]);
194 end
195 end
196 end
197 end