Mercurial > repos > public > sbplib

comparison +scheme/Laplace1d.m @ 1197:433c89bf19e0 feature/rv

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 07 Aug 2019 15:23:42 +0200
parents 33c378e508d2
children
comparison
equal deleted inserted replaced
1196:f6c571d8f22f 1197:433c89bf19e0
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','d','dirichlet'}
64 tuning = 1.1; 66 tuning = 1.1;
65 tau1 = -tuning/obj.gamm; 67 tau1 = -tuning/obj.gamm;
66 tau2 = 1; 68 tau2 = 1;
67 69
68 tau = tau1*e + tau2*d; 70 tau = tau1*e + tau2*d;
69 71
70 closure = obj.a*obj.Hi*tau*e'; 72 closure = obj.a*obj.Hi*tau*e';
71 penalty = obj.a*obj.Hi*tau; 73 penalty = -obj.a*obj.Hi*tau;
72 74
73 % Neumann boundary condition 75 % Neumann boundary condition
74 case {'N','neumann'} 76 case {'N','n','neumann'}
75 tau = -e; 77 tau = -e;
76 78
77 closure = obj.a*obj.Hi*tau*d'; 79 closure = obj.a*obj.Hi*tau*d';
78 penalty = -obj.a*obj.Hi*tau; 80 penalty = -obj.a*obj.Hi*tau;
79 81
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;
98 gamm_v = neighbour_scheme.gamm; 103 gamm_v = neighbour_scheme.gamm;
99 104
100 tuning = 1.1; 105 tuning = 1.1;
101 106
102 tau1 = -(a_u/gamm_u + a_v/gamm_v) * tuning; 107 tau1 = -1/4*(a_u/gamm_u + a_v/gamm_v) * tuning;
103 tau2 = 1/2*a_u; 108 tau2 = 1/2*a_u;
104 sig1 = -1/2; 109 sig1 = -1/2;
105 sig2 = 0; 110 sig2 = 0;
106 111
107 tau = tau1*e_u + tau2*d_u; 112 tau = tau1*e_u + tau2*d_u;
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();
133 end 155 end
134 156
135 end 157 end
136
137 methods(Static)
138 % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u
139 % and bound_v of scheme schm_v.
140 % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l')
141 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
142 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
143 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
144 end
145 end
146 end 158 end