Mercurial > repos > public > sbplib

comparison +scheme/bcSetup.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 1f6b2fb69225 6afd3dd3ed96
children 006defd0247b 57760d7088ad
comparison
equal deleted inserted replaced
796:aa1ed37a1b56 797:5cf9fdf4c98f
1 % function [closure, S] = bcSetup(diffOp, bc) 1 % function [closure, S] = bcSetup(diffOp, bc)
2 % Takes a diffOp and a cell array of boundary condition definitions. 2 % Takes a diffOp and a cell array of boundary condition definitions.
3 % Each bc is a struct with the fields 3 % Each bc is a struct with the fields
4 % * type -- Type of boundary condition 4 % * type -- Type of boundary condition
5 % * boundary -- Boundary identifier 5 % * boundary -- Boundary identifier
6 % * data -- A function_handle with time and space coordinates as a parameters, for example f(t,x,y) for a 2D problem 6 % * data -- A function_handle for a function which provides boundary data.(see below)
7 % Also takes S_sign which modifies the sign of S, [-1,1] 7 % Also takes S_sign which modifies the sign of S, [-1,1]
8 % Returns a closure matrix and a forcing function S 8 % Returns a closure matrix and a forcing function S.
9 function [closure, S] = bcSetup(diffOp, bc, S_sign) 9 %
10 % The boundary data function can either be a function of time or a function of time and space coordinates.
11 % In the case where it only depends on time it should return the data as grid function for the boundary.
12 % In the case where it also takes space coordinates the number of space coordinates should match the number of dimensions of the problem domain.
13 % For example in the 2D case: f(t,x,y).
14 function [closure, S] = bcSetup(diffOp, bcs, S_sign)
10 default_arg('S_sign', 1); 15 default_arg('S_sign', 1);
11 assertType(bc, 'cell'); 16 assertType(bcs, 'cell');
12 assert(S_sign == 1 || S_sign == -1, 'S_sign must be either 1 or -1'); 17 assert(S_sign == 1 || S_sign == -1, 'S_sign must be either 1 or -1');
13 18
19 verifyBcFormat(bcs, diffOp);
14 20
21 % Setup storage arrays
15 closure = spzeros(size(diffOp)); 22 closure = spzeros(size(diffOp));
16 penalties = {}; 23 gridData = {};
17 dataFunctions = {}; 24 symbolicData = {};
18 dataParams = {};
19 25
20 for i = 1:length(bc) 26 % Collect closures, penalties and data
21 assertType(bc{i}, 'struct'); 27 for i = 1:length(bcs)
22 [localClosure, penalty] = diffOp.boundary_condition(bc{i}.boundary, bc{i}.type); 28 [localClosure, penalty] = diffOp.boundary_condition(bcs{i}.boundary, bcs{i}.type);
23 closure = closure + localClosure; 29 closure = closure + localClosure;
24 30
25 if isempty(bc{i}.data) 31 [ok, isSym, data] = parseData(bcs{i}, penalty, diffOp.grid);
32
33 if ~ok
34 % There was no data
26 continue 35 continue
27 end 36 end
28 assertType(bc{i}.data, 'function_handle');
29 37
30 coord = diffOp.grid.getBoundary(bc{i}.boundary); 38 if isSym
31 assertNumberOfArguments(bc{i}.data, 1+size(coord,2)); 39 symbolicData{end+1} = data;
32 40 else
33 penalties{end+1} = penalty; 41 gridData{end+1} = data;
34 dataFunctions{end+1} = bc{i}.data; 42 end
35 dataParams{end+1} = num2cell(coord ,1);
36 end 43 end
37 44
45 % Setup penalty function
38 O = spzeros(size(diffOp),1); 46 O = spzeros(size(diffOp),1);
39 function v = S_fun(t) 47 function v = S_fun(t)
40 v = O; 48 v = O;
41 for i = 1:length(dataFunctions) 49 for i = 1:length(gridData)
42 v = v + penalties{i}*dataFunctions{i}(t, dataParams{i}{:}); 50 v = v + gridData{i}.penalty*gridData{i}.func(t);
51 end
52
53 for i = 1:length(symbolicData)
54 v = v + symbolicData{i}.penalty*symbolicData{i}.func(t, symbolicData{i}.coords{:});
43 end 55 end
44 56
45 v = S_sign * v; 57 v = S_sign * v;
46 end 58 end
47 S = @S_fun; 59 S = @S_fun;
48 end 60 end
61
62 function verifyBcFormat(bcs, diffOp)
63 for i = 1:length(bcs)
64 assertType(bcs{i}, 'struct');
65 assertStructFields(bcs{i}, {'type', 'boundary'});
66
67 if ~isfield(bcs{i}, 'data') || isempty(bcs{i}.data)
68 continue
69 end
70
71 if ~isa(bcs{i}.data, 'function_handle')
72 error('bcs{%d}.data should be a function of time or a function of time and space',i);
73 end
74
75 b = diffOp.grid.getBoundary(bcs{i}.boundary);
76
77 dim = size(b,2);
78
79 if nargin(bcs{i}.data) == 1
80 % Grid data (only function of time)
81 assertSize(bcs{i}.data(0), 1, size(b));
82 elseif nargin(bcs{i}.data) ~= 1+dim
83 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);
84 end
85 end
86 end
87
88 function [ok, isSym, dataStruct] = parseData(bc, penalty, grid)
89 if ~isfield(bc,'data') || isempty(bc.data)
90 ok = false;
91 return
92 end
93 ok = true;
94
95 nArg = nargin(bc.data);
96
97 if nArg > 1
98 % Symbolic data
99 isSym = true;
100 coord = grid.getBoundary(bc.boundary);
101 dataStruct.penalty = penalty;
102 dataStruct.func = bc.data;
103 dataStruct.coords = num2cell(coord, 1);
104 else
105 % Grid data
106 isSym = false;
107 dataStruct.penalty = penalty;
108 dataStruct.func = bcs{i}.data;
109 end
110 end