Mercurial > repos > public > sbplib

comparison +grid/Curvilinear.m @ 886:8894e9c49e40 feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default for latest changes
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 15 Nov 2018 16:36:21 -0800
parents 692bf61385c0
children
comparison
equal deleted inserted replaced
816:b5e5b195da1e 886:8894e9c49e40
1 classdef Curvilinear < grid.Structured & grid.Mapped
2 properties
3 logic % Grid of Logical domain
4 coords % N x D matrix with coordinates of each point in the physical domain
5 end
6
7 methods
8 % Creates a curvilinear grid.
9 % Ex: grid.Curvilinear(mapping, xi, eta, ...)
10 % mapping -- either a function handle, a matrix or a cell array with physical coordinates.
11 % A function handle should be a vector valued function of the coordinate mapping.
12 % A matrix should be a grid function (N*D x 1 vector) or a N x D
13 % A cell array should be a 1 x D cell array with either N x 1 vectors
14 % or matrices of the same dimesions as the logical grid.
15 % xi, eta, ... -- are the coordinate positions of the cartesian logical grid.
16 function obj = Curvilinear(mapping, varargin)
17 xi = varargin;
18 obj.logic = grid.Cartesian(xi{:});
19
20 % If mapping is a function evaluate it
21 if isa(mapping, 'function_handle')
22 if nargin(mapping) ~= length(varargin)
23 error('The dimension of the mapping does not match the dimension of the logical coordinates')
24 end
25 mapping = grid.evalOn(obj.logic, mapping);
26 end
27
28 D = obj.logic.D();
29 N = obj.logic.N();
30
31 obj.coords = zeros(N,D);
32
33 if iscell(mapping)
34 obj.coords = cellMappingToCoords(mapping, N, D, obj.logic.m);
35 elseif isnumeric(mapping)
36 obj.coords = matrixMappingToCoords(mapping, N, D);
37 else
38 error('grid:Curvilinear:Curvilinear','mapping must be a function handle, a matrix or a cell array.');
39 end
40 end
41
42 function m = size(obj)
43 m = obj.logic.size();
44 end
45
46 % logicalGrid returns the domain grid of the mapping.
47 function g = logicalGrid(obj)
48 g = obj.logic;
49 end
50
51 % mapping returns the mapped coordinates as a grid.Function
52 function m = mapping(obj);
53 m = obj.coords;
54 end
55
56 % n returns the number of points in the grid
57 function o = N(obj)
58 o = obj.logic.N();
59 end
60
61 % d returns the spatial dimension of the grid
62 function o = D(obj)
63 o = obj.logic.D();
64 end
65
66 % points returns a n x d matrix containing the coordinates for all points.
67 function X = points(obj)
68 X = obj.coords;
69 end
70
71 % Restricts the grid function gf on obj to the subgrid g.
72 function gf = restrictFunc(obj, gf, g)
73 gf = obj.logic.restrictFunc(gf, g.logic);
74 end
75
76 % Projects the grid function gf on obj to the grid g.
77 function gf = projectFunc(obj, gf, g)
78 gf = obj.logic.projectFunc(gf,g.logic);
79 end
80
81 function h = scaling(obj)
82 if isempty(obj.logic.h)
83 error('grid:Curvilinear:NoScalingSet','No scaling set');
84 end
85 h = obj.logic.h;
86 end
87
88 % Return the names of all boundaries in this grid.
89 function bs = getBoundaryNames(obj)
90 bs = obj.logic.getBoundaryNames();
91 end
92
93 % Return coordinates for the given boundary
94 function X = getBoundary(obj, name)
95 % In what dimension is the boundary?
96 switch name
97 case {'l', 'r', 'w', 'e'}
98 D = 1;
99 case {'s', 'n'}
100 D = 2;
101 case {'d', 'u'}
102 D = 3;
103 otherwise
104 error('not implemented');
105 end
106
107 % At what index is the boundary?
108 switch name
109 case {'l', 'w', 's', 'd'}
110 index = 1;
111 case {'r', 'e', 'n', 'u'}
112 index = obj.logic.m(D);
113 otherwise
114 error('not implemented');
115 end
116
117
118
119 I = cell(1, obj.D);
120 for i = 1:obj.D
121 if i == D
122 I{i} = index;
123 else
124 I{i} = ':';
125 end
126 end
127
128 % Calculate size of result:
129 m = obj.logic.m;
130 m(D) = [];
131 N = prod(m);
132
133 X = zeros(N, obj.D);
134
135 p = obj.points;
136 for i = 1:obj.D()
137 coordMat{i} = reshapeRowMaj(p(:,i), obj.logic.m);
138 end
139
140 for i = 1:length(coordMat)
141 Xtemp = coordMat{i}(I{:});
142 X(:,i) = reshapeRowMaj(Xtemp, [N,1]);
143 end
144 end
145 end
146 end
147
148
149 function coords = cellMappingToCoords(mapping, N, D, m)
150 if ~isequal(size(mapping),[1 D])
151 error('grid:Curvilinear:Curvilinear','The cell array must be a 1xD array.');
152 end
153
154 if isequal(size(mapping{1}),[N 1])
155 coords = cell2mat(mapping);
156 elseif isequal(size(mapping{1}), m)
157 for i = 1:length(mapping)
158 coords(:,i) = reshapeRowMaj(mapping{i}, [N 1]);
159 end
160 else
161 error('grid:Curvilinear:Curvilinear','The matrix must have size [N 1] or the same dimension as the grid. Actual: %s', toString(m));
162 end
163 end
164
165 function coords = matrixMappingToCoords(mapping, N, D)
166 if isequal(size(mapping), [N, D])
167 coords = mapping;
168 elseif isequal(size(mapping), [N*D, 1])
169 coords = reshapeRowMaj(mapping,[N D]);
170 else
171 error('grid:Curvilinear:Curvilinear','A matrix mapping must be of size [N D] or [N*D 1].');
172 end
173 end