Mercurial > repos > public > sbplib

comparison +multiblock/DiffOp.m @ 797:5cf9fdf4c98f feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with feature/grids and bugfix bcSetup
author Martin Almquist <malmquist@stanford.edu>
date Thu, 26 Jul 2018 10:53:05 -0700
parents aa1ed37a1b56 e05465aa2e25
children 57760d7088ad 7d4f57725192
comparison
equal deleted inserted replaced
796:aa1ed37a1b56 797:5cf9fdf4c98f
8 8
9 blockmatrixDiv 9 blockmatrixDiv
10 end 10 end
11 11
12 methods 12 methods
13 function obj = DiffOp(doHand, grid, order, doParam) 13 function obj = DiffOp(doHand, g, 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.
19 % grid -- a multiblock grid 19 % g -- a multiblock grid
20 % order -- integer specifying the order of accuracy 20 % order -- integer specifying the order of accuracy
21 % doParam -- may either be a cell array or a cell array of cell arrays 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 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 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, g, doParam);
30 30
31 nBlocks = grid.nBlocks(); 31 nBlocks = g.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);
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(g.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);
63 D{i,i} = obj.diffOps{i}.D; 63 D{i,i} = obj.diffOps{i}.D;
64 end 64 end
65 65
66 for i = 1:nBlocks 66 for i = 1:nBlocks
67 for j = 1:nBlocks 67 for j = 1:nBlocks
68 intf = grid.connections{i,j}; 68 intf = g.connections{i,j};
69 if isempty(intf) 69 if isempty(intf)
70 continue 70 continue
71 end 71 end
72 72
73 73
79 D{j,j} = D{j,j} + jj; 79 D{j,j} = D{j,j} + jj;
80 D{j,i} = D{j,i} + ji; 80 D{j,i} = D{j,i} + ji;
81 end 81 end
82 end 82 end
83 obj.D = blockmatrix.toMatrix(D); 83 obj.D = blockmatrix.toMatrix(D);
84 obj.grid = grid; 84 obj.grid = g;
85 85
86 86
87 function [getHand, getParam] = parseInput(doHand, grid, doParam) 87 function [getHand, getParam] = parseInput(doHand, g, doParam)
88 if ~isa(grid, 'multiblock.Grid') 88 if ~isa(g, 'multiblock.Grid')
89 error('multiblock:DiffOp:DiffOp:InvalidGrid', 'Requires a multiblock grid.'); 89 error('multiblock:DiffOp:DiffOp:InvalidGrid', 'Requires a multiblock grid.');
90 end 90 end
91 91
92 if iscell(doHand) && length(doHand) == grid.nBlocks() 92 if iscell(doHand) && length(doHand) == g.nBlocks()
93 getHand = @(i)doHand{i}; 93 getHand = @(i)doHand{i};
94 elseif isa(doHand, 'function_handle') 94 elseif isa(doHand, 'function_handle')
95 getHand = @(i)doHand; 95 getHand = @(i)doHand;
96 else 96 else
97 error('multiblock:DiffOp:DiffOp:InvalidGridDoHand', 'doHand must be a function handle or a cell array of length grid.nBlocks'); 97 error('multiblock:DiffOp:DiffOp:InvalidGridDoHand', 'doHand must be a function handle or a cell array of length grid.nBlocks');
107 return 107 return
108 end 108 end
109 109
110 % doParam is a non-empty cell-array 110 % doParam is a non-empty cell-array
111 111
112 if length(doParam) == grid.nBlocks() && all(cellfun(@iscell, doParam)) 112 if length(doParam) == g.nBlocks() && all(cellfun(@iscell, doParam))
113 % doParam is a cell-array of cell-arrays 113 % doParam is a cell-array of cell-arrays
114 getParam = @(i)doParam{i}; 114 getParam = @(i)doParam{i};
115 return 115 return
116 end 116 end
117 117
141 case 'multiblock.BoundaryGroup' 141 case 'multiblock.BoundaryGroup'
142 op = sparse(size(obj.D,1),0); 142 op = sparse(size(obj.D,1),0);
143 for i = 1:length(boundary) 143 for i = 1:length(boundary)
144 op = [op, obj.getBoundaryOperator(opName, boundary{i})]; 144 op = [op, obj.getBoundaryOperator(opName, boundary{i})];
145 end 145 end
146 otherwise
147 error('Unknown boundary indentifier')
148 end
149 end
150
151 function op = getBoundaryQuadrature(obj, boundary)
152 opName = 'H';
153 switch class(boundary)
154 case 'cell'
155 localOpName = [opName '_' boundary{2}];
156 blockId = boundary{1};
157 op = obj.diffOps{blockId}.(localOpName);
158
159 return
160 case 'multiblock.BoundaryGroup'
161 N = length(boundary);
162 H_bm = cell(N,N);
163 for i = 1:N
164 H_bm{i,i} = obj.getBoundaryQuadrature(boundary{i});
165 end
166 op = blockmatrix.toMatrix(H_bm);
146 otherwise 167 otherwise
147 error('Unknown boundary indentifier') 168 error('Unknown boundary indentifier')
148 end 169 end
149 end 170 end
150 171