Mercurial > repos > public > sbplib
changeset 186:1fc2eeb4f4e6 feature/grids
Moved multiblock grid to the multiblock package. Continued implementation. Added multiblock diffOp and removed some functions no longer needed.
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Thu, 03 Mar 2016 20:03:06 +0100 |
| parents | fad5e81389c1 |
| children | 770da87a6dc4 |
| files | +grid/Multiblock.m +multiblock/DiffOp.m +multiblock/DiffOpTest.m +multiblock/Grid.m +multiblock/GridTest.m +multiblock/gridVector1d.m +multiblock/gridVector2d.m +multiblock/solutionVector2cell.m |
| diffstat | 8 files changed, 243 insertions(+), 125 deletions(-) [+] |
line wrap: on
line diff
--- a/+grid/Multiblock.m Thu Mar 03 20:01:09 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,92 +0,0 @@ -classdef Multiblock < grid.Grid - properties - grids - connections - - nPoints - end - - % General multiblock grid - methods - - % grids -- cell array of N grids - % connections -- NxN cell matrix. connections{i,j} specifies the connection - % between block i and j. If it's empty there is no - % connection otherwise it's a 2-cell-vector with strings - % naming the boundaries to be connected. (inverted coupling?) - function obj = Multiblock(grids, connections) - obj.grids = grids; - obj.connections = connections; - - obj.nPoints = 0; - for i = 1:length(grids) - obj.nPoints = obj.nPoints + grids{i}.N(); - end - end - - function n = size(obj) - n = length(obj.grids); - end - - % n returns the number of points in the grid - function o = N(obj) - o = obj.nPoints; - end - - % d returns the spatial dimension of the grid - function o = D(obj) - o = obj.grids{1}.D(); - end - - % points returns a n x d matrix containing the coordinates for all points. - function X = points(obj) - X = []; - for i = 1:length(obj.grids) - X = [X; obj.grids{i}.points]; - end - end - - % Split a grid function on obj to a cell array of grid function on each block - function gfs = splitFunc(obj, gf) - nComponents = length(gf)/obj.nPoints; - nBlocks = length(obj.grids); - - % Collect number of points in each block - N = cell(1,nBlocks); - for i = 1:nBlocks - N{i} = obj.grids{i}.N(); - end - - gfs = mat2cell(gf, N, 1); - end - - % Restricts the grid function gf on obj to the subgrid g. - function gf = restrictFunc(obj, gf, g) - gfs = obj.splitFunc(gf); - - for i = 1:length(obj.grids) - gfs{i} = obj.grids{i}.restrictFunc(gfs{i}, g.grids{i}); - end - - gf = cell2mat(gfs); - end - - % Projects the grid function gf on obj to the grid g. - function o = projectFunc(obj, gf, g) - error('not implemented') - - p = g.points(); - o = zeros(length(p),1); - for i = 1:length(p) - I = whatGrid(p(i)); - o(i) = obj.grids{I}.projectFunc(gf, p(i)); - end - - - function I = whatGrid(p) - % Find what grid a point lies on - end - - end - end -end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/DiffOp.m Thu Mar 03 20:03:06 2016 +0100 @@ -0,0 +1,121 @@ +classdef DiffOp < scheme.Scheme + properties + grid + order + diffOps + D + H + end + + methods + function obj = DiffOp(doHand, grid, order, doParam) + % doHand -- may either be a function handle or a cell array of + % function handles for each grid. The function handle(s) + % should be on the form do = doHand(grid, order, ...) + % Additional parameters for each doHand may be provided in + % the doParam input. + % grid -- a multiblock grid + % order -- integer specifying the order of accuracy + % doParam -- may either be a cell array or a cell array of cell arrays + % for each block. If it is a cell array with length equal + % to the number of blocks then each element is sent to the + % corresponding function handle as extra parameters: + % doHand(..., doParam{i}{:}) Otherwise doParam is sent as + % extra parameters to all doHand: doHand(..., doParam{:}) + default_arg('doParam', []) + + [getHand, getParam] = parseInput(doHand, grid, doParam); + + nBlocks = grid.nBlocks(); + + obj.order = order; + + % Create the diffOps for each block + obj.diffOps = cell(1, nBlocks); + for i = 1:nBlocks + h = getHand(i); + p = getParam(i); + obj.diffOps{i} = h(grid.grid{i}, order, p{:}); + end + + + % Build the norm matrix + H = cell(nBlocks, nBlocks); + for i = 1:nBlocks + H{i,i} = obj.diffOps{i}.H; + end + obj.H = cell2sparse(H); + + + % Build the differentiation matrix + D = cell(nBlocks, nBlocks); + for i = 1:nBlocks + D{i,i} = obj.diffOps{i}.D; + end + + for i = 1:nBlocks + for j = i:nBlocks + intf = grid.connections{i,j}; + if isempty(intf) + continue + end + + [ii, ij] = obj.diffOps{i}.inteface_coupling(intf{1}, obj.diffOps{j}, intf{2}); + D{i,i} = D{i,i} + ii; + D{i,j} = D{i,j} + ij; + + [jj, ji] = obj.diffOps{j}.inteface_coupling(intf{2}, obj.diffOps{i}, intf{1}); + D{j,j} = D{j,j} + jj; + D{j,i} = D{j,i} + ji; + end + end + obj.D = cell2sparse(D); + + % end + + function [getHand, getParam] = parseInput(doHand, grid, doParam) + if ~isa(grid, 'multiblock.Grid') + error('multiblock:DiffOp:DiffOp:InvalidGrid', 'Requires a multiblock grid.'); + end + + if iscell(doHand) && length(doHand) == grid.nBlocks() + getHand = @(i)doHand{i}; + elseif isa(doHand, 'function_handle') + getHand = @(i)doHand; + else + error('multiblock:DiffOp:DiffOp:InvalidGridDoHand', 'doHand must be a function handle or a cell array of length grid.nBlocks'); + end + + if isempty(doParam) + getParam = @(i){}; + elseif iscell(doParam) && length(doParam) == grid.nBlocks() + getParam = @(i)doParam{i}; + else + getParam = @(i)doParam; + end + end + end + + function ops = splitOp(obj, op) + % Splits a matrix operator into a cell-matrix of matrix operators for + % each grid. + ops = sparse2cell(op, obj.NNN); + end + + function m = boundary_condition(obj,boundary,type,data) + + end + + function m = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + + end + + + function N = size(obj) + N = 0; + for i = 1:length(diffOps) + N = N + diffOps.size(); + end + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/DiffOpTest.m Thu Mar 03 20:03:06 2016 +0100 @@ -0,0 +1,18 @@ +function tests = DiffOpTest() + tests = functiontests(localfunctions); +end + +function testCreation(testCase) + g = multiblock.Grid({},{}); + doHand = @(grid,order)[]; + order = 0; + do = multiblock.DiffOp(doHand, g, order); +end + + + + +% function do = mockDiffOp() +% do.H = 1; +% do.D = 1; +% end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/Grid.m Thu Mar 03 20:03:06 2016 +0100 @@ -0,0 +1,97 @@ +classdef Grid < grid.Grid + properties + grids + connections + + nPoints + end + + % General multiblock grid + methods + + % grids -- cell array of N grids + % connections -- NxN upper triangular cell matrix. connections{i,j} + % specifies the connection between block i and j. If + % it's empty there is no connection otherwise it's a 2 + % -cell-vector with strings naming the boundaries to be + % connected. (inverted coupling?) + function obj = Grid(grids, connections) + obj.grids = grids; + obj.connections = connections; + + obj.nPoints = 0; + for i = 1:length(grids) + obj.nPoints = obj.nPoints + grids{i}.N(); + end + end + + function n = size(obj) + n = length(obj.grids); + end + + % n returns the number of points in the grid + function o = N(obj) + o = obj.nPoints; + end + + function n = nBlocks(obj) + n = length(obj.grids); + end + + % d returns the spatial dimension of the grid + function o = D(obj) + o = obj.grids{1}.D(); + end + + % points returns a n x d matrix containing the coordinates for all points. + function X = points(obj) + X = []; + for i = 1:length(obj.grids) + X = [X; obj.grids{i}.points]; + end + end + + % Split a grid function on obj to a cell array of grid function on each block + function gfs = splitFunc(obj, gf) + nComponents = length(gf)/obj.nPoints; + nBlocks = length(obj.grids); + + % Collect number of points in each block + N = cell(1,nBlocks); + for i = 1:nBlocks + N{i} = obj.grids{i}.N(); + end + + gfs = mat2cell(gf, N, 1); + end + + % Restricts the grid function gf on obj to the subgrid g. + function gf = restrictFunc(obj, gf, g) + gfs = obj.splitFunc(gf); + + for i = 1:length(obj.grids) + gfs{i} = obj.grids{i}.restrictFunc(gfs{i}, g.grids{i}); + end + + gf = cell2mat(gfs); + end + + % Projects the grid function gf on obj to the grid g. + function o = projectFunc(obj, gf, g) + error('not implemented') + + p = g.points(); + o = zeros(length(p),1); + for i = 1:length(p) + I = whatGrid(p(i)); + o(i) = obj.grids{I}.projectFunc(gf, p(i)); + end + + + function I = whatGrid(p) + % Find what grid a point lies on + end + + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/GridTest.m Thu Mar 03 20:03:06 2016 +0100 @@ -0,0 +1,7 @@ +function tests = GridTest() + tests = functiontests(localfunctions); +end + +function testCreation(testCase) + g = multiblock.Grid({},{}); +end \ No newline at end of file
--- a/+multiblock/gridVector1d.m Thu Mar 03 20:01:09 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,8 +0,0 @@ -function x = gridVector1d(schms) - n = length(schms); - x = cell(n,1); - - for i = 1:n - x{i} = schms{i}.x; - end -end \ No newline at end of file
--- a/+multiblock/gridVector2d.m Thu Mar 03 20:01:09 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,10 +0,0 @@ -function [x,y] = gridVector2d(schms) - n = length(schms); - x = cell(n,1); - y = cell(n,1); - - for i = 1:n - x{i} = schms{i}.x; - y{i} = schms{i}.y; - end -end \ No newline at end of file
--- a/+multiblock/solutionVector2cell.m Thu Mar 03 20:01:09 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,15 +0,0 @@ -function v_cell = solutionVector2cell(schms,v,components) - if nargin == 2 - components = 1; - end - - N = length(schms); - v_cell = cell(N,1); - - i_start = 1; - for i =1:N - i_end = i_start + schms{i}.size()*components - 1; - v_cell{i} = v(i_start:i_end); - i_start = i_end+1; - end -end \ No newline at end of file