Mercurial > repos > public > sbplib

comparison +grid/Curvilinear.m @ 427:a613960a157b feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merged with feature/beams
author Ylva Rydin <ylva.rydin@telia.com>
date Thu, 26 Jan 2017 15:59:25 +0100
parents 4f935415700e
children 692bf61385c0
comparison
equal deleted inserted replaced
426:29944ea7674b 427:a613960a157b
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 mapping = grid.evalOn(obj.logic, mapping);
23 end
24
25 D = obj.logic.D();
26 N = obj.logic.N();
27
28 obj.coords = zeros(N,D);
29
30 if iscell(mapping)
31 obj.coords = cellMappingToCoords(mapping, N, D, obj.logic.m);
32 elseif isnumeric(mapping)
33 obj.coords = matrixMappingToCoords(mapping, N, D);
34 else
35 error('grid:Curvilinear:Curvilinear','mapping must be a function handle, a matrix or a cell array.');
36 end
37 end
38
39 function m = size(obj)
40 m = obj.logic.size();
41 end
42
43 % logicalGrid returns the domain grid of the mapping.
44 function g = logicalGrid(obj)
45 g = obj.logic;
46 end
47
48 % mapping returns the mapped coordinates as a grid.Function
49 function m = mapping(obj);
50 m = obj.coords;
51 end
52
53 % n returns the number of points in the grid
54 function o = N(obj)
55 o = obj.logic.N();
56 end
57
58 % d returns the spatial dimension of the grid
59 function o = D(obj)
60 o = obj.logic.D();
61 end
62
63 % points returns a n x d matrix containing the coordinates for all points.
64 function X = points(obj)
65 X = obj.coords;
66 end
67
68 % Restricts the grid function gf on obj to the subgrid g.
69 function gf = restrictFunc(obj, gf, g)
70 gf = obj.logic.restrictFunc(gf, g.logic);
71 end
72
73 % Projects the grid function gf on obj to the grid g.
74 function gf = projectFunc(obj, gf, g)
75 gf = obj.logic.projectFunc(gf,g.logic);
76 end
77
78 function h = scaling(obj)
79 if isempty(obj.logic.h)
80 error('grid:Curvilinear:NoScalingSet','No scaling set');
81 end
82 h = obj.logic.h;
83 end
84
85 % Return the names of all boundaries in this grid.
86 function bs = getBoundaryNames(obj)
87 bs = obj.logic.getBoundaryNames();
88 end
89
90 % Return coordinates for the given boundary
91 function X = getBoundary(obj, name)
92 % In what dimension is the boundary?
93 switch name
94 case {'l', 'r', 'w', 'e'}
95 D = 1;
96 case {'s', 'n'}
97 D = 2;
98 case {'d', 'u'}
99 D = 3;
100 otherwise
101 error('not implemented');
102 end
103
104 % At what index is the boundary?
105 switch name
106 case {'l', 'w', 's', 'd'}
107 index = 1;
108 case {'r', 'e', 'n', 'u'}
109 index = obj.logic.m(D);
110 otherwise
111 error('not implemented');
112 end
113
114
115
116 I = cell(1, obj.D);
117 for i = 1:obj.D
118 if i == D
119 I{i} = index;
120 else
121 I{i} = ':';
122 end
123 end
124
125 % Calculate size of result:
126 m = obj.logic.m;
127 m(D) = [];
128 N = prod(m);
129
130 X = zeros(N, obj.D);
131
132 p = obj.points;
133 for i = 1:obj.D()
134 coordMat{i} = reshapeRowMaj(p(:,i), obj.logic.m);
135 end
136
137 for i = 1:length(coordMat)
138 Xtemp = coordMat{i}(I{:});
139 X(:,i) = reshapeRowMaj(Xtemp, [N,1]);
140 end
141 end
142 end
143 end
144
145
146 function coords = cellMappingToCoords(mapping, N, D, m)
147 if ~isequal(size(mapping),[1 D])
148 error('grid:Curvilinear:Curvilinear','The cell array must be a 1xD array.');
149 end
150
151 if isequal(size(mapping{1}),[N 1])
152 coords = cell2mat(mapping);
153 elseif isequal(size(mapping{1}), m)
154 for i = 1:length(mapping)
155 coords(:,i) = reshapeRowMaj(mapping{i}, [N 1]);
156 end
157 else
158 error('grid:Curvilinear:Curvilinear','The matrix must have size [N 1] or the same dimension as the grid. Actual: %s', toString(m));
159 end
160 end
161
162 function coords = matrixMappingToCoords(mapping, N, D)
163 if isequal(size(mapping), [N, D])
164 coords = mapping;
165 elseif isequal(size(mapping), [N*D, 1])
166 coords = reshapeRowMaj(mapping,[N D]);
167 else
168 error('grid:Curvilinear:Curvilinear','A matrix mapping must be of size [N D] or [N*D 1].');
169 end
170 end