Mercurial > repos > public > sbplib

comparison +scheme/Utux.m @ 1072:6468a5f6ec79 feature/grids/LaplaceSquared

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 12 Feb 2019 17:12:42 +0100
parents 0c504a21432d
children 433c89bf19e0
comparison
equal deleted inserted replaced
1071:92cb03e64ca4 1072:6468a5f6ec79
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
14 e_r 14 e_r
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, opSet)
21 default_arg('a',1); 21 default_arg('opSet',@sbp.D2Standard);
22
23 %Old operators
24 % [x, h] = util.get_grid(xlim{:},m);
25 %ops = sbp.Ordinary(m,h,order);
26
27
28 switch operator
29 case 'NonEquidistant'
30 ops = sbp.D1Nonequidistant(m,xlim,order);
31 obj.D1 = ops.D1;
32 case 'Standard'
33 ops = sbp.D2Standard(m,xlim,order);
34 obj.D1 = ops.D1;
35 case 'Upwind'
36 ops = sbp.D1Upwind(m,xlim,order);
37 obj.D1 = ops.Dm;
38 otherwise
39 error('Unvalid operator')
40 end
41 obj.x=ops.x;
42 22
43 23 m = g.size();
24 xl = g.getBoundary('l');
25 xr = g.getBoundary('r');
26 xlim = {xl, xr};
27
28 ops = opSet(m, xlim, order);
29 obj.D1 = ops.D1;
30
31 obj.grid = g;
32
44 obj.H = ops.H; 33 obj.H = ops.H;
45 obj.Hi = ops.HI; 34 obj.Hi = ops.HI;
46 35
47 obj.e_l = ops.e_l; 36 obj.e_l = ops.e_l;
48 obj.e_r = ops.e_r; 37 obj.e_r = ops.e_r;
49 obj.D=obj.D1; 38 obj.D = -obj.D1;
50 39
51 obj.m = m; 40 obj.m = m;
52 obj.h = ops.h; 41 obj.h = ops.h;
53 obj.order = order; 42 obj.order = order;
54 obj.x = ops.x;
55 43
56 end 44 end
57 % Closure functions return the opertors applied to the own doamin to close the boundary 45 % 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. 46 % 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'. 47 % 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. 48 % 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. 49 % 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. 50 % neighbour_scheme is an instance of Scheme that should be interfaced to.
63 % neighbour_boundary is a string specifying which boundary to interface to. 51 % neighbour_boundary is a string specifying which boundary to interface to.
64 function [closure, penalty] = boundary_condition(obj,boundary,type,data) 52 function [closure, penalty] = boundary_condition(obj,boundary,type)
65 default_arg('type','neumann'); 53 default_arg('type','dirichlet');
66 default_arg('data',0); 54 tau =-1*obj.e_l;
67 tau =-1*obj.e_l; 55 closure = obj.Hi*tau*obj.e_l';
68 closure = obj.Hi*tau*obj.e_l'; 56 penalty = -obj.Hi*tau;
69 penalty = 0*obj.e_l; 57
70
71 end 58 end
72 59
73 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 60 function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type)
74 error('An interface function does not exist yet'); 61 switch boundary
62 % Upwind coupling
63 case {'l','left'}
64 tau = -1*obj.e_l;
65 closure = obj.Hi*tau*obj.e_l';
66 penalty = -obj.Hi*tau*neighbour_scheme.e_r';
67 case {'r','right'}
68 tau = 0*obj.e_r;
69 closure = obj.Hi*tau*obj.e_r';
70 penalty = -obj.Hi*tau*neighbour_scheme.e_l';
71 end
72
75 end 73 end
76 74
75 % Returns the boundary operator op for the boundary specified by the string boundary.
76 % op -- string
77 % boundary -- string
78 function o = getBoundaryOperator(obj, op, boundary)
79 assertIsMember(op, {'e'})
80 assertIsMember(boundary, {'l', 'r'})
81
82 o = obj.([op, '_', boundary]);
83 end
84
85 % Returns square boundary quadrature matrix, of dimension
86 % corresponding to the number of boundary points
87 %
88 % boundary -- string
89 % Note: for 1d diffOps, the boundary quadrature is the scalar 1.
90 function H_b = getBoundaryQuadrature(obj, boundary)
91 assertIsMember(boundary, {'l', 'r'})
92
93 H_b = 1;
94 end
95
77 function N = size(obj) 96 function N = size(obj)
78 N = obj.m; 97 N = obj.m;
79 end 98 end
80 99
81 end 100 end
82
83 methods(Static)
84 % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u
85 % and bound_v of scheme schm_v.
86 % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l')
87 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);
89 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
90 end
91 end
92 end 101 end