Mercurial > repos > public > sbplib

comparison +scheme/Laplace1d.m @ 1108:5ec23b9bf360 feature/laplace_curvilinear_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Martin Almquist <malmquist@stanford.edu>
date Wed, 10 Apr 2019 11:00:27 -0700
parents a22b23098021
children 33c378e508d2
comparison
equal deleted inserted replaced
1087:867307f4d80f 1108:5ec23b9bf360
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] = obj.getBoundaryOperator({'e', 'd'}, boundary); 59 e = obj.getBoundaryOperator('e', boundary);
60 d = obj.getBoundaryOperator('d', boundary);
60 s = obj.getBoundarySign(boundary); 61 s = obj.getBoundarySign(boundary);
61 62
62 switch type 63 switch type
63 % Dirichlet boundary condition 64 % Dirichlet boundary condition
64 case {'D','dirichlet'} 65 case {'D','d','dirichlet'}
65 tuning = 1.1; 66 tuning = 1.1;
66 tau1 = -tuning/obj.gamm; 67 tau1 = -tuning/obj.gamm;
67 tau2 = 1; 68 tau2 = 1;
68 69
69 tau = tau1*e + tau2*d; 70 tau = tau1*e + tau2*d;
70 71
71 closure = obj.a*obj.Hi*tau*e'; 72 closure = obj.a*obj.Hi*tau*e';
72 penalty = obj.a*obj.Hi*tau; 73 penalty = obj.a*obj.Hi*tau;
73 74
74 % Neumann boundary condition 75 % Neumann boundary condition
75 case {'N','neumann'} 76 case {'N','n','neumann'}
76 tau = -e; 77 tau = -e;
77 78
78 closure = obj.a*obj.Hi*tau*d'; 79 closure = obj.a*obj.Hi*tau*d';
79 penalty = -obj.a*obj.Hi*tau; 80 penalty = -obj.a*obj.Hi*tau;
80 81
85 end 86 end
86 87
87 function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type) 88 function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type)
88 % u denotes the solution in the own domain 89 % u denotes the solution in the own domain
89 % v denotes the solution in the neighbour domain 90 % v denotes the solution in the neighbour domain
90 [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary); 91 e_u = obj.getBoundaryOperator('e', boundary);
92 d_u = obj.getBoundaryOperator('d', boundary);
91 s_u = obj.getBoundarySign(boundary); 93 s_u = obj.getBoundarySign(boundary);
92 94
93 [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary); 95 e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
96 d_v = neighbour_scheme.getBoundaryOperator('d', neighbour_boundary);
94 s_v = neighbour_scheme.getBoundarySign(neighbour_boundary); 97 s_v = neighbour_scheme.getBoundarySign(neighbour_boundary);
95 98
96 a_u = obj.a; 99 a_u = obj.a;
97 a_v = neighbour_scheme.a; 100 a_v = neighbour_scheme.a;
98 101
99 gamm_u = obj.gamm; 102 gamm_u = obj.gamm;
100 gamm_v = neighbour_scheme.gamm; 103 gamm_v = neighbour_scheme.gamm;
101 104
102 tuning = 1.1; 105 tuning = 1.1;
103 106
104 tau1 = -(a_u/gamm_u + a_v/gamm_v) * tuning; 107 tau1 = -1/4*(a_u/gamm_u + a_v/gamm_v) * tuning;
105 tau2 = 1/2*a_u; 108 tau2 = 1/2*a_u;
106 sig1 = -1/2; 109 sig1 = -1/2;
107 sig2 = 0; 110 sig2 = 0;
108 111
109 tau = tau1*e_u + tau2*d_u; 112 tau = tau1*e_u + tau2*d_u;
112 closure = obj.Hi*( tau*e_u' + sig*a_u*d_u'); 115 closure = obj.Hi*( tau*e_u' + sig*a_u*d_u');
113 penalty = obj.Hi*(-tau*e_v' + sig*a_v*d_v'); 116 penalty = obj.Hi*(-tau*e_v' + sig*a_v*d_v');
114 end 117 end
115 118
116 % Returns the boundary operator op for the boundary specified by the string boundary. 119 % Returns the boundary operator op for the boundary specified by the string boundary.
117 % op -- string or a cell array of strings 120 % op -- string
118 % boundary -- string 121 % boundary -- string
119 function varargout = getBoundaryOperator(obj, op, boundary) 122 function o = getBoundaryOperator(obj, op, boundary)
123 assertIsMember(op, {'e', 'd'})
124 assertIsMember(boundary, {'l', 'r'})
120 125
121 if ~iscell(op) 126 o = obj.([op, '_', boundary]);
122 op = {op}; 127 end
123 end
124 128
125 for i = 1:numel(op) 129 % Returns square boundary quadrature matrix, of dimension
126 switch op{i} 130 % corresponding to the number of boundary points
127 case 'e' 131 %
128 switch boundary 132 % boundary -- string
129 case 'l' 133 % Note: for 1d diffOps, the boundary quadrature is the scalar 1.
130 e = obj.e_l; 134 function H_b = getBoundaryQuadrature(obj, boundary)
131 case 'r' 135 assertIsMember(boundary, {'l', 'r'})
132 e = obj.e_r;
133 otherwise
134 error('No such boundary: boundary = %s',boundary);
135 end
136 varargout{i} = e;
137 136
138 case 'd' 137 H_b = 1;
139 switch boundary
140 case 'l'
141 d = obj.d_l;
142 case 'r'
143 d = obj.d_r;
144 otherwise
145 error('No such boundary: boundary = %s',boundary);
146 end
147 varargout{i} = d;
148 end
149 end
150 end 138 end
151 139
152 % Returns the boundary sign. The right boundary is considered the positive boundary 140 % Returns the boundary sign. The right boundary is considered the positive boundary
153 % boundary -- string 141 % boundary -- string
154 function s = getBoundarySign(obj, boundary) 142 function s = getBoundarySign(obj, boundary)
143 assertIsMember(boundary, {'l', 'r'})
144
155 switch boundary 145 switch boundary
156 case {'r'} 146 case {'r'}
157 s = 1; 147 s = 1;
158 case {'l'} 148 case {'l'}
159 s = -1; 149 s = -1;
160 otherwise
161 error('No such boundary: boundary = %s',boundary);
162 end 150 end
163 end 151 end
164 152
165 function N = size(obj) 153 function N = size(obj)
166 N = obj.grid.size(); 154 N = obj.grid.size();
167 end 155 end
168 156
169 end 157 end
170
171 methods(Static)
172 % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u
173 % and bound_v of scheme schm_v.
174 % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l')
175 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
176 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
177 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
178 end
179 end
180 end 158 end