Mercurial > repos > public > sbplib

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