Mercurial > repos > public > sbplib

comparison +multiblock/DiffOp.m @ 707:0de70ec8bf60 feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merge with feature/optim
author Ylva Rydin <ylva.rydin@telia.com>
date Fri, 10 Nov 2017 14:22:56 +0100
parents 111fcbcff2e9
children
comparison
equal deleted inserted replaced
696:7c16b5af8d98 707:0de70ec8bf60
3 grid 3 grid
4 order 4 order
5 diffOps 5 diffOps
6 D 6 D
7 H 7 H
8 8
9 blockmatrixDiv 9 blockmatrixDiv
10 end 10 end
11 11
12 methods 12 methods
13 function obj = DiffOp(doHand, grid, order, doParam,timeDep) 13 function obj = DiffOp(doHand, grid, order, doParam)
14 % doHand -- may either be a function handle or a cell array of 14 % doHand -- may either be a function handle or a cell array of
15 % function handles for each grid. The function handle(s) 15 % function handles for each grid. The function handle(s)
16 % should be on the form do = doHand(grid, order, ...) 16 % should be on the form do = doHand(grid, order, ...)
17 % Additional parameters for each doHand may be provided in 17 % Additional parameters for each doHand may be provided in
18 % the doParam input. 18 % the doParam input.
23 % to the number of blocks then each element is sent to the 23 % to the number of blocks then each element is sent to the
24 % corresponding function handle as extra parameters: 24 % corresponding function handle as extra parameters:
25 % doHand(..., doParam{i}{:}) Otherwise doParam is sent as 25 % doHand(..., doParam{i}{:}) Otherwise doParam is sent as
26 % extra parameters to all doHand: doHand(..., doParam{:}) 26 % extra parameters to all doHand: doHand(..., doParam{:})
27 default_arg('doParam', []) 27 default_arg('doParam', [])
28 28
29 [getHand, getParam] = parseInput(doHand, grid, doParam); 29 [getHand, getParam] = parseInput(doHand, grid, doParam);
30 30
31 nBlocks = grid.nBlocks(); 31 nBlocks = grid.nBlocks();
32 32
33 obj.order = order; 33 obj.order = order;
34 34
35 % Create the diffOps for each block 35 % Create the diffOps for each block
36 obj.diffOps = cell(1, nBlocks); 36 obj.diffOps = cell(1, nBlocks);
37 for i = 1:nBlocks 37 for i = 1:nBlocks
38 h = getHand(i); 38 h = getHand(i);
39 p = getParam(i); 39 p = getParam(i);
40 if ~iscell(p) 40 if ~iscell(p)
41 p = {p}; 41 p = {p};
42 end 42 end
43 obj.diffOps{i} = h(grid.grids{i}, order, p{:}); 43 obj.diffOps{i} = h(grid.grids{i}, order, p{:});
44 end 44 end
45 45
46 46
47 % Build the norm matrix 47 % Build the norm matrix
48 H = cell(nBlocks, nBlocks); 48 H = cell(nBlocks, nBlocks);
49 for i = 1:nBlocks 49 for i = 1:nBlocks
50 H{i,i} = obj.diffOps{i}.H; 50 H{i,i} = obj.diffOps{i}.H;
51 end 51 end
52 obj.H = blockmatrix.toMatrix(H); 52 obj.H = blockmatrix.toMatrix(H);
53 53
54 54
55 % Build the differentiation matrix 55 % Build the differentiation matrix
56 switch timeDep 56 obj.blockmatrixDiv = {grid.Ns, grid.Ns};
57 case {'n','no','N','No'} 57 D = blockmatrix.zero(obj.blockmatrixDiv);
58 obj.blockmatrixDiv = {grid.Ns, grid.Ns}; 58 for i = 1:nBlocks
59 D = blockmatrix.zero(obj.blockmatrixDiv); 59 D{i,i} = obj.diffOps{i}.D;
60 for i = 1:nBlocks 60 end
61 D{i,i} = obj.diffOps{i}.D; 61
62 for i = 1:nBlocks
63 for j = 1:nBlocks
64 intf = grid.connections{i,j};
65 if isempty(intf)
66 continue
62 end 67 end
63 68
64 for i = 1:nBlocks 69
65 for j = 1:nBlocks 70 [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2});
66 intf = grid.connections{i,j}; 71 D{i,i} = D{i,i} + ii;
67 if isempty(intf) 72 D{i,j} = D{i,j} + ij;
68 continue 73
69 end 74 [jj, ji] = obj.diffOps{j}.interface(intf{2}, obj.diffOps{i}, intf{1});
70 75 D{j,j} = D{j,j} + jj;
71 76 D{j,i} = D{j,i} + ji;
72 [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2}); 77 end
73 D{i,i} = D{i,i} + ii; 78 end
74 D{i,j} = D{i,j} + ij; 79 obj.D = blockmatrix.toMatrix(D);
75 80
76 [jj, ji] = obj.diffOps{j}.interface(intf{2}, obj.diffOps{i}, intf{1}); 81
77 D{j,j} = D{j,j} + jj;
78 D{j,i} = D{j,i} + ji;
79 end
80 end
81 obj.D = blockmatrix.toMatrix(D);
82 case {'y','yes','Y','Yes'}
83 for i = 1:nBlocks
84 D{i,i} = @(t)obj.diffOps{i}.D(t);
85 end
86
87 for i = 1:nBlocks
88 for j = 1:nBlocks
89 intf = grid.connections{i,j};
90 if isempty(intf)
91 continue
92 end
93
94 [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2});
95 D{i,i} = @(t)D{i,i}(t) + ii(t);
96 D{i,j} = ij;
97
98 [jj, ji] = obj.diffOps{j}.interface(intf{2}, obj.diffOps{i}, intf{1});
99 D{j,j} = @(t)D{j,j}(t) + jj(t);
100 D{j,i} = ji;
101 end
102 end
103 obj.D = D;
104 end
105
106
107 function [getHand, getParam] = parseInput(doHand, grid, doParam) 82 function [getHand, getParam] = parseInput(doHand, grid, doParam)
108 if ~isa(grid, 'multiblock.Grid') 83 if ~isa(grid, 'multiblock.Grid')
109 error('multiblock:DiffOp:DiffOp:InvalidGrid', 'Requires a multiblock grid.'); 84 error('multiblock:DiffOp:DiffOp:InvalidGrid', 'Requires a multiblock grid.');
110 end 85 end
111 86
112 if iscell(doHand) && length(doHand) == grid.nBlocks() 87 if iscell(doHand) && length(doHand) == grid.nBlocks()
113 getHand = @(i)doHand{i}; 88 getHand = @(i)doHand{i};
114 elseif isa(doHand, 'function_handle') 89 elseif isa(doHand, 'function_handle')
115 getHand = @(i)doHand; 90 getHand = @(i)doHand;
116 else 91 else
117 error('multiblock:DiffOp:DiffOp:InvalidGridDoHand', 'doHand must be a function handle or a cell array of length grid.nBlocks'); 92 error('multiblock:DiffOp:DiffOp:InvalidGridDoHand', 'doHand must be a function handle or a cell array of length grid.nBlocks');
118 end 93 end
119 94
120 if isempty(doParam) 95 if isempty(doParam)
121 getParam = @(i){}; 96 getParam = @(i){};
122 return 97 return
123 end 98 end
124 99
125 if ~iscell(doParam) 100 if ~iscell(doParam)
126 getParam = @(i)doParam; 101 getParam = @(i)doParam;
127 return 102 return
128 end 103 end
129 104
130 % doParam is a non-empty cell-array 105 % doParam is a non-empty cell-array
131 106
132 if length(doParam) == grid.nBlocks() && all(cellfun(@iscell, doParam)) 107 if length(doParam) == grid.nBlocks() && all(cellfun(@iscell, doParam))
133 % doParam is a cell-array of cell-arrays 108 % doParam is a cell-array of cell-arrays
134 getParam = @(i)doParam{i}; 109 getParam = @(i)doParam{i};
135 return 110 return
136 end 111 end
137 112
138 getParam = @(i)doParam; 113 getParam = @(i)doParam;
139 end 114 end
140 end 115 end
141 116
142 function ops = splitOp(obj, op) 117 function ops = splitOp(obj, op)
143 % Splits a matrix operator into a cell-matrix of matrix operators for 118 % Splits a matrix operator into a cell-matrix of matrix operators for
144 % each grid. 119 % each grid.
145 ops = sparse2cell(op, obj.NNN); 120 ops = sparse2cell(op, obj.NNN);
146 end 121 end
147 122
148 function op = getBoundaryOperator(obj, op, boundary) 123 % Get a boundary operator specified by opName for the given boundary/BoundaryGroup
149 if iscell(boundary) 124 function op = getBoundaryOperator(obj, opName, boundary)
150 localOpName = [op '_' boundary{2}]; 125 switch class(boundary)
151 blockId = boundary{1}; 126 case 'cell'
152 localOp = obj.diffOps{blockId}.(localOpName); 127 localOpName = [opName '_' boundary{2}];
153 128 blockId = boundary{1};
154 div = {obj.blockmatrixDiv{1}, size(localOp,2)}; 129 localOp = obj.diffOps{blockId}.(localOpName);
155 blockOp = blockmatrix.zero(div); 130
156 blockOp{blockId,1} = localOp; 131 div = {obj.blockmatrixDiv{1}, size(localOp,2)};
157 op = blockmatrix.toMatrix(blockOp); 132 blockOp = blockmatrix.zero(div);
158 return 133 blockOp{blockId,1} = localOp;
159 else 134 op = blockmatrix.toMatrix(blockOp);
160 % Boundary är en sträng med en boundary group i. 135 return
161 end 136 case 'multiblock.BoundaryGroup'
162 end 137 op = sparse(size(obj.D,1),0);
163 138 for i = 1:length(boundary)
139 op = [op, obj.getBoundaryOperator(opName, boundary{i})];
140 end
141 otherwise
142 error('Unknown boundary indentifier')
143 end
144 end
145
164 % Creates the closure and penalty matrix for a given boundary condition, 146 % Creates the closure and penalty matrix for a given boundary condition,
165 % boundary -- the name of the boundary on the form {id,name} where 147 % boundary -- the name of the boundary on the form {id,name} where
166 % id is the number of a block and name is the name of a 148 % id is the number of a block and name is the name of a
167 % boundary of that block example: {1,'s'} or {3,'w'} 149 % boundary of that block example: {1,'s'} or {3,'w'}. It
150 % can also be a boundary group
168 function [closure, penalty] = boundary_condition(obj, boundary, type) 151 function [closure, penalty] = boundary_condition(obj, boundary, type)
152 switch class(boundary)
153 case 'cell'
154 [closure, penalty] = obj.singleBoundaryCondition(boundary, type);
155 case 'multiblock.BoundaryGroup'
156 [n,m] = size(obj.D);
157 closure = sparse(n,m);
158 penalty = sparse(n,0);
159 for i = 1:length(boundary)
160 [closurePart, penaltyPart] = obj.boundary_condition(boundary{i}, type);
161 closure = closure + closurePart;
162 penalty = [penalty, penaltyPart];
163 end
164 otherwise
165 error('Unknown boundary indentifier')
166 end
167
168 end
169
170 function [closure, penalty] = singleBoundaryCondition(obj, boundary, type)
169 I = boundary{1}; 171 I = boundary{1};
170 name = boundary{2}; 172 name = boundary{2};
171 173
172 % Get the closure and penaly matrices 174 % Get the closure and penaly matrices
173 [blockClosure, blockPenalty] = obj.diffOps{I}.boundary_condition(name, type); 175 [blockClosure, blockPenalty] = obj.diffOps{I}.boundary_condition(name, type);
174 176
175 % Expand to matrix for full domain. 177 % Expand to matrix for full domain.
176 div = obj.blockmatrixDiv; 178 div = obj.blockmatrixDiv;
177 if ~iscell(blockClosure) 179 if ~iscell(blockClosure)
178 temp = blockmatrix.zero(div); 180 temp = blockmatrix.zero(div);
179 temp{I,I} = blockClosure; 181 temp{I,I} = blockClosure;
183 temp = blockmatrix.zero(div); 185 temp = blockmatrix.zero(div);
184 temp{I,I} = blockClosure{i}; 186 temp{I,I} = blockClosure{i};
185 closure{i} = blockmatrix.toMatrix(temp); 187 closure{i} = blockmatrix.toMatrix(temp);
186 end 188 end
187 end 189 end
188 190
189 div{2} = size(blockPenalty, 2); % Penalty is a column vector
190 if ~iscell(blockPenalty) 191 if ~iscell(blockPenalty)
192 div{2} = size(blockPenalty, 2); % Penalty is a column vector
191 p = blockmatrix.zero(div); 193 p = blockmatrix.zero(div);
192 p{I} = blockPenalty; 194 p{I} = blockPenalty;
193 penalty = blockmatrix.toMatrix(p); 195 penalty = blockmatrix.toMatrix(p);
194 else 196 else
197 % TODO: used by beam equation, should be eliminated. SHould only set one BC per call
195 for i = 1:length(blockPenalty) 198 for i = 1:length(blockPenalty)
199 div{2} = size(blockPenalty{i}, 2); % Penalty is a column vector
196 p = blockmatrix.zero(div); 200 p = blockmatrix.zero(div);
197 p{I} = blockPenalty{i}; 201 p{I} = blockPenalty{i};
198 penalty{i} = blockmatrix.toMatrix(p); 202 penalty{i} = blockmatrix.toMatrix(p);
199 end 203 end
200 end 204 end
201 end 205 end
202 206
203 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 207 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
204 208 error('not implemented')
205 end 209 end
206 210
207 % Size returns the number of degrees of freedom 211 % Size returns the number of degrees of freedom
208 function N = size(obj) 212 function N = size(obj)
209 N = 0; 213 N = 0;
210 for i = 1:length(obj.diffOps) 214 for i = 1:length(obj.diffOps)
211 N = N + obj.diffOps{i}.size(); 215 N = N + obj.diffOps{i}.size();