Mercurial > repos > public > sbplib

comparison +scheme/Laplace1d.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 ae4b090b5299
comparison
equal deleted inserted replaced
1071:92cb03e64ca4 1072:6468a5f6ec79
1 classdef Laplace1d < scheme.Scheme
2 properties
3 grid
4 order % Order accuracy for the approximation
5
6 D % non-stabalized scheme operator
7 H % Discrete norm
8 M % Derivative norm
9 a
10
11 D2
12 Hi
13 e_l
14 e_r
15 d_l
16 d_r
17 gamm
18 end
19
20 methods
21 function obj = Laplace1d(grid, order, a)
22 default_arg('a', 1);
23
24 assertType(grid, 'grid.Cartesian');
25
26 ops = sbp.D2Standard(grid.size(), grid.lim{1}, order);
27
28 obj.D2 = sparse(ops.D2);
29 obj.H = sparse(ops.H);
30 obj.Hi = sparse(ops.HI);
31 obj.M = sparse(ops.M);
32 obj.e_l = sparse(ops.e_l);
33 obj.e_r = sparse(ops.e_r);
34 obj.d_l = -sparse(ops.d1_l);
35 obj.d_r = sparse(ops.d1_r);
36
37
38 obj.grid = grid;
39 obj.order = order;
40
41 obj.a = a;
42 obj.D = a*obj.D2;
43
44 obj.gamm = grid.h*ops.borrowing.M.S;
45 end
46
47
48 % Closure functions return the opertors applied to the own doamin to close the boundary
49 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
50 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
51 % type is a string specifying the type of boundary condition if there are several.
52 % data is a function returning the data that should be applied at the boundary.
53 % neighbour_scheme is an instance of Scheme that should be interfaced to.
54 % neighbour_boundary is a string specifying which boundary to interface to.
55 function [closure, penalty] = boundary_condition(obj,boundary,type,data)
56 default_arg('type','neumann');
57 default_arg('data',0);
58
59 e = obj.getBoundaryOperator('e', boundary);
60 d = obj.getBoundaryOperator('d', boundary);
61 s = obj.getBoundarySign(boundary);
62
63 switch type
64 % Dirichlet boundary condition
65 case {'D','dirichlet'}
66 tuning = 1.1;
67 tau1 = -tuning/obj.gamm;
68 tau2 = 1;
69
70 tau = tau1*e + tau2*d;
71
72 closure = obj.a*obj.Hi*tau*e';
73 penalty = obj.a*obj.Hi*tau;
74
75 % Neumann boundary condition
76 case {'N','neumann'}
77 tau = -e;
78
79 closure = obj.a*obj.Hi*tau*d';
80 penalty = -obj.a*obj.Hi*tau;
81
82 % Unknown, boundary condition
83 otherwise
84 error('No such boundary condition: type = %s',type);
85 end
86 end
87
88 function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type)
89 % u denotes the solution in the own 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);
94
95 e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
96 d_v = neighbour_scheme.getBoundaryOperator('d', neighbour_boundary);
97 s_v = neighbour_scheme.getBoundarySign(neighbour_boundary);
98
99 a_u = obj.a;
100 a_v = neighbour_scheme.a;
101
102 gamm_u = obj.gamm;
103 gamm_v = neighbour_scheme.gamm;
104
105 tuning = 1.1;
106
107 tau1 = -(a_u/gamm_u + a_v/gamm_v) * tuning;
108 tau2 = 1/2*a_u;
109 sig1 = -1/2;
110 sig2 = 0;
111
112 tau = tau1*e_u + tau2*d_u;
113 sig = sig1*e_u + sig2*d_u;
114
115 closure = obj.Hi*( tau*e_u' + sig*a_u*d_u');
116 penalty = obj.Hi*(-tau*e_v' + sig*a_v*d_v');
117 end
118
119 % Returns the boundary operator op for the boundary specified by the string boundary.
120 % op -- string
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
145 switch boundary
146 case {'r'}
147 s = 1;
148 case {'l'}
149 s = -1;
150 end
151 end
152
153 function N = size(obj)
154 N = obj.grid.size();
155 end
156
157 end
158 end