|
705
|
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
|
