Mercurial > repos > public > sbplib
changeset 791:d6100c4fe3b6 feature/grids
Merge with bcSetupExperiment
author | Jonatan Werpers <jonatan@werpers.com> |
---|---|
date | Wed, 25 Jul 2018 09:38:06 -0700 |
parents | 085fc0fe537f (current diff) 6afd3dd3ed96 (diff) |
children | c7c622e26a53 77aa40da2008 |
files | |
diffstat | 1 files changed, 68 insertions(+), 33 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/bcSetup.m Tue Jul 24 20:17:25 2018 -0700 +++ b/+scheme/bcSetup.m Wed Jul 25 09:38:06 2018 -0700 @@ -11,49 +11,34 @@ % In the case where it only depends on time it should return the data as grid function for the boundary. % In the case where it also takes space coordinates the number of space coordinates should match the number of dimensions of the problem domain. % For example in the 2D case: f(t,x,y). -function [closure, S] = bcSetup(diffOp, bc, S_sign) +function [closure, S] = bcSetup(diffOp, bcs, S_sign) default_arg('S_sign', 1); - assertType(bc, 'cell'); + assertType(bcs, 'cell'); assert(S_sign == 1 || S_sign == -1, 'S_sign must be either 1 or -1'); + verifyBcFormat(bcs, diffOp); % Setup storage arrays closure = spzeros(size(diffOp)); - gridDataPenalties = {}; - gridDataFunctions = {}; - symbolicDataPenalties = {}; - symbolicDataFunctions = {}; - symbolicDataCoords = {}; + gridData = {}; + symbolicData = {}; % Collect closures, penalties and data - for i = 1:length(bc) - assertType(bc{i}, 'struct'); - [localClosure, penalty] = diffOp.boundary_condition(bc{i}.boundary, bc{i}.type); + for i = 1:length(bcs) + [localClosure, penalty] = diffOp.boundary_condition(bcs{i}.boundary, bcs{i}.type); closure = closure + localClosure; - if ~isfield(bc{i},'data') || isempty(bc{i}.data) - % Skip to next loop if there is no data + [ok, isSym, data] = parseData(bcs{i}, penalty, diffOp.grid) + + if ~ok + % There was no data continue end - assertType(bc{i}.data, 'function_handle'); - % Find dimension - dim = size(diffOp.grid.getBoundary(bc{i}.boundary), 2); - - if nargin(bc{i}.data) == 1 - % Grid data - boundarySize = [size(diffOp.grid.getBoundary(bc{i}.boundary),1),1]; - assertSize(bc{i}.data(0), boundarySize); % Eval for t = 0 and make sure the function returns a grid vector of the correct size. - gridDataPenalties{end+1} = penalty; - gridDataFunctions{end+1} = bc{i}.data; - elseif nargin(bc{i}.data) == 1+dim - % Symbolic data - coord = diffOp.grid.getBoundary(bc{i}.boundary); - symbolicDataPenalties{end+1} = penalty; - symbolicDataFunctions{end+1} = bc{i}.data; - symbolicDataCoords{end+1} = num2cell(coord ,1); + if isSym + gridData{end+1} = data; else - error('sbplib:scheme:bcSetup:DataWrongNumberOfArguments', 'bc{%d}.data has the wrong number of input arguments. Must be either only time or time and space.', i); + symbolicData{end+1} = data; end end @@ -61,15 +46,65 @@ O = spzeros(size(diffOp),1); function v = S_fun(t) v = O; - for i = 1:length(gridDataFunctions) - v = v + gridDataPenalties{i}*gridDataFunctions{i}(t); + for i = 1:length(gridData) + v = v + gridData{i}.penalty*gridData{i}.func(t); end - for i = 1:length(symbolicDataFunctions) - v = v + symbolicDataPenalties{i}*symbolicDataFunctions{i}(t, symbolicDataCoords{i}{:}); + for i = 1:length(symbolicData) + v = v + symbolicData{i}.penalty*symbolicData{i}.func(t, symbolicData{i}.coords{:}); end v = S_sign * v; end S = @S_fun; end + +function verifyBcFormat(bcs, diffOp) + for i = 1:length(bcs) + assertType(bcs{i}, 'struct'); + assertStructFields(bcs{i}, {'type', 'boundary'}); + + if ~isfield(bcs{i}, 'data') || isempty(bcs{i}.data) + continue + end + + if ~isa(bcs{i}.data, 'function_handle') + error('bcs{%d}.data should be a function of time or a function of time and space',i); + end + + b = diffOp.grid.getBoundary(bc.boundart); + + dim = size(b,2); + + if nargin(bc.data) == 1 + % Grid data (only function of time) + assertSize(bc.data(0), 1, size(b)); + elseif nargin(bc.data) ~= 1+dim + error('sbplib:scheme:bcSetup:DataWrongNumberOfArguments', 'bcs{%d}.data has the wrong number of input arguments. Must be either only time or time and space.', i); + end + end +end + +function [ok, isSym, dataStruct] = parseData(bc, penalty, grid) + if ~isfield(bc,'data') || isempty(bc.data) + ok = false; + return + end + ok = true; + + nArg = nargin(bc.data); + + if nArg > 1 + % Symbolic data + isSym = true; + coord = grid.getBoundary(bc.boundary); + dataStruct.penalty = penalty; + dataStruct.func = bc.data; + dataStruct.coords = num2cell(coord, 1); + else + % Grid data + isSym = false; + dataStruct.penalty = penalty; + dataStruct.func = bcs{i}.data; + end +end