Mercurial > repos > public > sbplib

comparison +scheme/Laplace1d.m @ 1054:77676c26056d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merged in feature/getBoundaryOp (pull request #12) Feature/getBoundaryOp
author Jonatan Werpers <jonatan.werpers@it.uu.se>
date Thu, 24 Jan 2019 08:29:50 +0000
parents 0c504a21432d
children ae4b090b5299
comparison
equal deleted inserted replaced
1043:c12b84fe9b00 1054:77676c26056d
54 % neighbour_boundary is a string specifying which boundary to interface to. 54 % neighbour_boundary is a string specifying which boundary to interface to.
55 function [closure, penalty] = boundary_condition(obj,boundary,type,data) 55 function [closure, penalty] = boundary_condition(obj,boundary,type,data)
56 default_arg('type','neumann'); 56 default_arg('type','neumann');
57 default_arg('data',0); 57 default_arg('data',0);
58 58
59 [e,d,s] = obj.get_boundary_ops(boundary); 59 e = obj.getBoundaryOperator('e', boundary);
60 d = obj.getBoundaryOperator('d', boundary);
61 s = obj.getBoundarySign(boundary);
60 62
61 switch type 63 switch type
62 % Dirichlet boundary condition 64 % Dirichlet boundary condition
63 case {'D','dirichlet'} 65 case {'D','dirichlet'}
64 tuning = 1.1; 66 tuning = 1.1;
84 end 86 end
85 87
86 function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type) 88 function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type)
87 % u denotes the solution in the own domain 89 % u denotes the solution in the own domain
88 % v denotes the solution in the neighbour domain 90 % v denotes the solution in the neighbour domain
91 e_u = obj.getBoundaryOperator('e', boundary);
92 d_u = obj.getBoundaryOperator('d', boundary);
93 s_u = obj.getBoundarySign(boundary);
89 94
90 [e_u,d_u,s_u] = obj.get_boundary_ops(boundary); 95 e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
91 [e_v,d_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 96 d_v = neighbour_scheme.getBoundaryOperator('d', neighbour_boundary);
92 97 s_v = neighbour_scheme.getBoundarySign(neighbour_boundary);
93 98
94 a_u = obj.a; 99 a_u = obj.a;
95 a_v = neighbour_scheme.a; 100 a_v = neighbour_scheme.a;
96 101
97 gamm_u = obj.gamm; 102 gamm_u = obj.gamm;
109 114
110 closure = obj.Hi*( tau*e_u' + sig*a_u*d_u'); 115 closure = obj.Hi*( tau*e_u' + sig*a_u*d_u');
111 penalty = obj.Hi*(-tau*e_v' + sig*a_v*d_v'); 116 penalty = obj.Hi*(-tau*e_v' + sig*a_v*d_v');
112 end 117 end
113 118
114 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. 119 % Returns the boundary operator op for the boundary specified by the string boundary.
115 % The right boundary is considered the positive boundary 120 % op -- string
116 function [e,d,s] = get_boundary_ops(obj,boundary) 121 % boundary -- string
122 function o = getBoundaryOperator(obj, op, boundary)
123 assertIsMember(op, {'e', 'd'})
124 assertIsMember(boundary, {'l', 'r'})
125
126 o = obj.([op, '_', boundary])
127 end
128
129 % Returns square boundary quadrature matrix, of dimension
130 % corresponding to the number of boundary points
131 %
132 % boundary -- string
133 % Note: for 1d diffOps, the boundary quadrature is the scalar 1.
134 function H_b = getBoundaryQuadrature(obj, boundary)
135 assertIsMember(boundary, {'l', 'r'})
136
137 H_b = 1;
138 end
139
140 % Returns the boundary sign. The right boundary is considered the positive boundary
141 % boundary -- string
142 function s = getBoundarySign(obj, boundary)
143 assertIsMember(boundary, {'l', 'r'})
144
117 switch boundary 145 switch boundary
118 case 'l' 146 case {'r'}
119 e = obj.e_l; 147 s = 1;
120 d = obj.d_l; 148 case {'l'}
121 s = -1; 149 s = -1;
122 case 'r'
123 e = obj.e_r;
124 d = obj.d_r;
125 s = 1;
126 otherwise
127 error('No such boundary: boundary = %s',boundary);
128 end 150 end
129 end 151 end
130 152
131 function N = size(obj) 153 function N = size(obj)
132 N = obj.grid.size(); 154 N = obj.grid.size();