Mercurial > repos > public > sbplib

comparison +scheme/Burgers1d.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 f6c571d8f22f
children 68ee061639a1
comparison
equal deleted inserted replaced
1196:f6c571d8f22f 1197:433c89bf19e0
66 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other domain. 66 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other domain.
67 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 67 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
68 % type is a string specifying the type of boundary condition if there are several. 68 % type is a string specifying the type of boundary condition if there are several.
69 function [closure, penalty] = boundary_condition(obj, boundary, type) 69 function [closure, penalty] = boundary_condition(obj, boundary, type)
70 default_arg('type','dirichlet'); 70 default_arg('type','dirichlet');
71 [e, index, s] = obj.get_boundary_ops(boundary); 71 s = obj.getBoundarySign(boundary);
72 e = obj.getBoundaryOperator('e', boundary);
73 index = obj.getBoundaryIndex(boundary);
72 switch type 74 switch type
73 % Stable dirchlet-like boundary conditions (u+-abs(u))*u/3 75 % Stable dirchlet-like boundary conditions (u+-abs(u))*u/3
74 % with +- at left/right boundaries 76 % with +- at left/right boundaries
75 case {'D', 'd', 'dirichlet', 'Dirichlet'} 77 case {'D', 'd', 'dirichlet', 'Dirichlet'}
76 % tau = s*e; 78 % tau = s*e;
77 % closure = @(v) obj.Hi*tau*(((v(index)-s*abs(v(index)))/3)*v(index)); 79 % closure = @(v) obj.Hi*tau*(((v(index)-s*abs(v(index)))/3)*v(index));
78 % penalty = -obj.Hi*tau; 80 % penalty = -obj.Hi*tau;
79 81
80 magnitude = 2/3; 82 penalty_parameter = 1/3;
81 tau = @(v) s*magnitude*obj.Hi*e*(v(index)-s*abs(v(index)))/2; 83 tau = @(v) s*penalty_parameter*obj.Hi*e*(v(index)-s*abs(v(index)))/2;
82 closure = @(v) tau(v)*v(index); 84 closure = @(v) tau(v)*v(index);
83 penalty = @(v) -tau(v); 85 penalty = @(v) -tau(v);
84 otherwise 86 otherwise
85 error('No such boundary condition: type = %s',type); 87 error('No such boundary condition: type = %s',type);
86 end 88 end
87 end 89 end
88 90
89 % Returns the boundary ops, boundary index and sign for the boundary specified by the string boundary. 91
90 % The right boundary is considered the positive boundary 92 % Returns the boundary sign. The right boundary is considered the positive boundary
91 function [e, index, s] = get_boundary_ops(obj,boundary) 93 % boundary -- string
94 function s = getBoundarySign(obj, boundary)
95 assertIsMember(boundary, {'l', 'r'})
96
92 switch boundary 97 switch boundary
93 case {'l','L','left','Left'} 98 case {'r'}
94 e = obj.e_l; 99 s = 1;
100 case {'l'}
101 s = -1;
102 end
103 end
104
105 % Returns the boundary operator op for the boundary specified by the string boundary.
106 % op -- string
107 % boundary -- string
108 function o = getBoundaryOperator(obj, op, boundary)
109 assertIsMember(op, {'e'})
110 assertIsMember(boundary, {'l', 'r'})
111
112 o = obj.([op, '_', boundary]);
113 end
114
115 % Returns square boundary quadrature matrix, of dimension
116 % corresponding to the number of boundary points
117 %
118 % boundary -- string
119 % Note: for 1d diffOps, the boundary quadrature is the scalar 1.
120 function H_b = getBoundaryQuadrature(obj, boundary)
121 assertIsMember(boundary, {'l', 'r'})
122 H_b = 1;
123 end
124
125 % Returns the boundary index. The right boundary has the last index
126 % boundary -- string
127 function index = getBoundaryIndex(obj, boundary)
128 assertIsMember(boundary, {'l', 'r'})
129 switch boundary
130 case {'r'}
131 index = length(obj.e_r);
132 case {'l'}
95 index = 1; 133 index = 1;
96 s = -1;
97 case {'r','R','right','Right'}
98 e = obj.e_r;
99 index = length(e);
100 s = 1;
101 otherwise
102 error('No such boundary: boundary = %s',boundary);
103 end 134 end
104 end 135 end
105 136
106 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 137 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
107 error('An interface function does not exist yet'); 138 error('An interface function does not exist yet');