Mercurial > repos > public > sbplib
annotate +multiblock/DiffOpTimeDep.m @ 713:348d5bcf7daf feature/quantumTriangles
Merge with feature/frids
| author | Ylva Rydin <ylva.rydin@telia.com> |
|---|---|
| date | Tue, 20 Feb 2018 15:00:30 +0100 |
| parents | 397d1b22cc37 |
| children |
| rev | line source |
|---|---|
| 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(); | |
|
711
397d1b22cc37
Changed closure and penalty ao that sparse matrices are created in DiffOpTimeDep
Ylva Rydin <ylva.rydin@telia.com>
parents:
708
diff
changeset
|
166 [n,m] = size(obj.D{1,1}(0)); |
| 705 | 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 | |
|
711
397d1b22cc37
Changed closure and penalty ao that sparse matrices are created in DiffOpTimeDep
Ylva Rydin <ylva.rydin@telia.com>
parents:
708
diff
changeset
|
173 closure{j,i} = @(t)0*sparse(n,m); |
|
397d1b22cc37
Changed closure and penalty ao that sparse matrices are created in DiffOpTimeDep
Ylva Rydin <ylva.rydin@telia.com>
parents:
708
diff
changeset
|
174 penalty{j,i} = @(t)0*sparse(n,1); |
| 705 | 175 end |
| 176 end | |
| 177 | |
| 178 | |
| 179 for i = 1:length(boundary) | |
| 708 | 180 boundaryPart = boundary{i}; |
| 181 block = boundaryPart{1}; | |
| 182 [closurePart, penaltyPart] = obj.boundary_condition(boundaryPart, type); | |
| 183 closure{block,block} = @(t)closure{block,block}(t) + closurePart(t); | |
| 184 penalty{block,block} = @(t)penalty{block,block}(t) + penaltyPart(t); | |
| 705 | 185 end |
| 186 otherwise | |
| 187 error('Unknown boundary indentifier') | |
| 188 end | |
| 189 | |
| 190 end | |
| 191 | |
| 192 function [blockClosure, blockPenalty] = singleBoundaryCondition(obj, boundary, type) | |
| 193 I = boundary{1}; | |
| 194 name = boundary{2}; | |
| 195 % Get the closure and penaly matrices | |
| 196 [blockClosure, blockPenalty] = obj.diffOps{I}.boundary_condition(name, type); | |
| 197 | |
| 198 end | |
| 199 | |
| 200 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) | |
| 201 error('not implemented') | |
| 202 end | |
| 203 | |
| 204 % Size returns the number of degrees of freedom | |
| 205 function N = size(obj) | |
| 206 N = 0; | |
| 207 for i = 1:length(obj.diffOps) | |
| 208 N = N + obj.diffOps{i}.size(); | |
| 209 end | |
| 210 end | |
| 211 end | |
| 212 end |