Mercurial > repos > public > sbplib

comparison +scheme/Utux.m @ 756:f891758ad7a4 feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with feature/utux2d.
author Martin Almquist <malmquist@stanford.edu>
date Sat, 16 Jun 2018 14:30:45 -0700
parents efe2dbf9796e
children 459eeb99130f
comparison
equal deleted inserted replaced
755:14f0058356f2 756:f891758ad7a4
1 classdef Utux < scheme.Scheme 1 classdef Utux < scheme.Scheme
2 properties 2 properties
3 m % Number of points in each direction, possibly a vector 3 m % Number of points in each direction, possibly a vector
4 h % Grid spacing 4 h % Grid spacing
5 x % Grid 5 grid % Grid
6 order % Order accuracy for the approximation 6 order % Order accuracy for the approximation
7 7
8 H % Discrete norm 8 H % Discrete norm
9 D 9 D
10 10
15 v0 15 v0
16 end 16 end
17 17
18 18
19 methods 19 methods
20 function obj = Utux(m,xlim,order,operator) 20 function obj = Utux(g ,order, operator)
21 default_arg('a',1); 21 default_arg('operator','Standard');
22 22
23 %Old operators 23 m = g.size();
24 % [x, h] = util.get_grid(xlim{:},m); 24 xl = g.getBoundary('l');
25 %ops = sbp.Ordinary(m,h,order); 25 xr = g.getBoundary('r');
26 26 xlim = {xl, xr};
27 27
28 switch operator 28 switch operator
29 case 'NonEquidistant' 29 % case 'NonEquidistant'
30 ops = sbp.D1Nonequidistant(m,xlim,order); 30 % ops = sbp.D1Nonequidistant(m,xlim,order);
31 obj.D1 = ops.D1; 31 % obj.D1 = ops.D1;
32 case 'Standard' 32 case 'Standard'
33 ops = sbp.D2Standard(m,xlim,order); 33 ops = sbp.D2Standard(m,xlim,order);
34 obj.D1 = ops.D1; 34 obj.D1 = ops.D1;
35 case 'Upwind' 35 % case 'Upwind'
36 ops = sbp.D1Upwind(m,xlim,order); 36 % ops = sbp.D1Upwind(m,xlim,order);
37 obj.D1 = ops.Dm; 37 % obj.D1 = ops.Dm;
38 otherwise 38 otherwise
39 error('Unvalid operator') 39 error('Unvalid operator')
40 end 40 end
41 obj.x=ops.x; 41
42 obj.grid = g;
42 43
43
44 obj.H = ops.H; 44 obj.H = ops.H;
45 obj.Hi = ops.HI; 45 obj.Hi = ops.HI;
46 46
47 obj.e_l = ops.e_l; 47 obj.e_l = ops.e_l;
48 obj.e_r = ops.e_r; 48 obj.e_r = ops.e_r;
49 obj.D=obj.D1; 49 obj.D = -obj.D1;
50 50
51 obj.m = m; 51 obj.m = m;
52 obj.h = ops.h; 52 obj.h = ops.h;
53 obj.order = order; 53 obj.order = order;
54 obj.x = ops.x;
55 54
56 end 55 end
57 % Closure functions return the opertors applied to the own doamin to close the boundary 56 % Closure functions return the opertors applied to the own doamin to close the boundary
58 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. 57 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
59 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 58 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
60 % type is a string specifying the type of boundary condition if there are several. 59 % type is a string specifying the type of boundary condition if there are several.
61 % data is a function returning the data that should be applied at the boundary. 60 % data is a function returning the data that should be applied at the boundary.
62 % neighbour_scheme is an instance of Scheme that should be interfaced to. 61 % neighbour_scheme is an instance of Scheme that should be interfaced to.
63 % neighbour_boundary is a string specifying which boundary to interface to. 62 % neighbour_boundary is a string specifying which boundary to interface to.
64 function [closure, penalty] = boundary_condition(obj,boundary,type,data) 63 function [closure, penalty] = boundary_condition(obj,boundary,type)
65 default_arg('type','neumann'); 64 default_arg('type','dirichlet');
66 default_arg('data',0);
67 tau =-1*obj.e_l; 65 tau =-1*obj.e_l;
68 closure = obj.Hi*tau*obj.e_l'; 66 closure = obj.Hi*tau*obj.e_l';
69 penalty = 0*obj.e_l; 67 penalty = -obj.Hi*tau;
70 68
71 end 69 end
72 70
73 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 71 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
74 error('An interface function does not exist yet'); 72 switch boundary
73 % Upwind coupling
74 case {'l','left'}
75 tau = -1*obj.e_l;
76 closure = obj.Hi*tau*obj.e_l';
77 penalty = -obj.Hi*tau*neighbour_scheme.e_r';
78 case {'r','right'}
79 tau = 0*obj.e_r;
80 closure = obj.Hi*tau*obj.e_r';
81 penalty = -obj.Hi*tau*neighbour_scheme.e_l';
82 end
83
75 end 84 end
76 85
77 function N = size(obj) 86 function N = size(obj)
78 N = obj.m; 87 N = obj.m;
79 end 88 end
80 89
81 end 90 end
82 91
83 methods(Static) 92 methods(Static)
84 % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u 93 % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u
85 % and bound_v of scheme schm_v. 94 % and bound_v of scheme schm_v.
86 % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') 95 % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l')
87 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) 96 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
88 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); 97 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
89 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); 98 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
90 end 99 end
91 end 100 end