Mercurial > repos > public > sbplib

annotate +scheme/Hypsyst2d.m @ 1198:2924b3a9b921 feature/d2_compatible

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add OpSet for fully compatible D2Variable, created from regular D2Variable by replacing d1 by first row of D1. Formal reduction by one order of accuracy at the boundary point.
author Martin Almquist <malmquist@stanford.edu>
date Fri, 16 Aug 2019 14:30:28 -0700
parents 8d73fcdb07a5
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
295
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
1 classdef Hypsyst2d < scheme.Scheme
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
2 properties
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
3 m % Number of points in each direction, possibly a vector
293
2d604d16842c Works with varying coefficients and char boundary condition
Ylva Rydin <ylva.rydin@telia.com>
parents: 292
diff changeset
4 n %size of system
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
5 h % Grid spacing
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
6 x,y % Grid
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
7 X,Y % Values of x and y for each grid point
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
8 order % Order accuracy for the approximation
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
9
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
10 D % non-stabalized scheme operator
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
11 A, B, E %Coefficient matrices
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
12
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
13 H % Discrete norm
291
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
14 % Norms in the x and y directions
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
15 Hxi,Hyi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
16 I_x,I_y, I_N
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
17 e_w, e_e, e_s, e_n
297
cd30b22cee56 Have tried to make a curvelinear sheme for hypsysts. Does not really work yet...
Ylva Rydin <ylva.rydin@telia.com>
parents: 296
diff changeset
18 params %parameters for the coeficient matrice
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
19 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
20
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
21 methods
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
22 %Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Eu
295
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
23 function obj = Hypsyst2d(m, lim, order, A, B, E, params)
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
24 default_arg('E', [])
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
25 xlim = lim{1};
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
26 ylim = lim{2};
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
27
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
28 if length(m) == 1
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
29 m = [m m];
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
30 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
31
297
cd30b22cee56 Have tried to make a curvelinear sheme for hypsysts. Does not really work yet...
Ylva Rydin <ylva.rydin@telia.com>
parents: 296
diff changeset
32 obj.A=A;
cd30b22cee56 Have tried to make a curvelinear sheme for hypsysts. Does not really work yet...
Ylva Rydin <ylva.rydin@telia.com>
parents: 296
diff changeset
33 obj.B=B;
cd30b22cee56 Have tried to make a curvelinear sheme for hypsysts. Does not really work yet...
Ylva Rydin <ylva.rydin@telia.com>
parents: 296
diff changeset
34 obj.E=E;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
35
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
36 m_x = m(1);
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
37 m_y = m(2);
295
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
38 obj.params = params;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
39
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
40 ops_x = sbp.D2Standard(m_x,xlim,order);
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
41 ops_y = sbp.D2Standard(m_y,ylim,order);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
42
295
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
43 obj.x = ops_x.x;
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
44 obj.y = ops_y.x;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
45
291
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
46 obj.X = kr(obj.x,ones(m_y,1));
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
47 obj.Y = kr(ones(m_x,1),obj.y);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
48
297
cd30b22cee56 Have tried to make a curvelinear sheme for hypsysts. Does not really work yet...
Ylva Rydin <ylva.rydin@telia.com>
parents: 296
diff changeset
49 Aevaluated = obj.evaluateCoefficientMatrix(A, obj.X, obj.Y);
cd30b22cee56 Have tried to make a curvelinear sheme for hypsysts. Does not really work yet...
Ylva Rydin <ylva.rydin@telia.com>
parents: 296
diff changeset
50 Bevaluated = obj.evaluateCoefficientMatrix(B, obj.X, obj.Y);
cd30b22cee56 Have tried to make a curvelinear sheme for hypsysts. Does not really work yet...
Ylva Rydin <ylva.rydin@telia.com>
parents: 296
diff changeset
51 Eevaluated = obj.evaluateCoefficientMatrix(E, obj.X, obj.Y);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
52
297
cd30b22cee56 Have tried to make a curvelinear sheme for hypsysts. Does not really work yet...
Ylva Rydin <ylva.rydin@telia.com>
parents: 296
diff changeset
53 obj.n = length(A(obj.params,0,0));
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
54
295
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
55 I_n = eye(obj.n);I_x = speye(m_x);
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
56 obj.I_x = I_x;
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
57 I_y = speye(m_y);
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
58 obj.I_y = I_y;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
59
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
60
295
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
61 D1_x = kr(I_n, ops_x.D1, I_y);
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
62 obj.Hxi = kr(I_n, ops_x.HI, I_y);
297
cd30b22cee56 Have tried to make a curvelinear sheme for hypsysts. Does not really work yet...
Ylva Rydin <ylva.rydin@telia.com>
parents: 296
diff changeset
63 D1_y = kr(I_n, I_x, ops_y.D1);
cd30b22cee56 Have tried to make a curvelinear sheme for hypsysts. Does not really work yet...
Ylva Rydin <ylva.rydin@telia.com>
parents: 296
diff changeset
64 obj.Hyi = kr(I_n, I_x, ops_y.HI);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
65
295
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
66 obj.e_w = kr(I_n, ops_x.e_l, I_y);
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
67 obj.e_e = kr(I_n, ops_x.e_r, I_y);
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
68 obj.e_s = kr(I_n, I_x, ops_y.e_l);
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
69 obj.e_n = kr(I_n, I_x, ops_y.e_r);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
70
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
71 obj.m = m;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
72 obj.h = [ops_x.h ops_y.h];
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
73 obj.order = order;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
74
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
75 obj.D = -Aevaluated*D1_x-Bevaluated*D1_y-Eevaluated;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
76
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
77 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
78
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
79 % Closure functions return the opertors applied to the own doamin to close the boundary
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
80 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
81 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
82 % type is a string specifying the type of boundary condition if there are several.
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
83 % data is a function returning the data that should be applied at the boundary.
293
2d604d16842c Works with varying coefficients and char boundary condition
Ylva Rydin <ylva.rydin@telia.com>
parents: 292
diff changeset
84 function [closure, penalty] = boundary_condition(obj,boundary,type,L)
294
8ff6ec6249e8 "General" boundary conditions implemented
Ylva Rydin <ylva.rydin@telia.com>
parents: 293
diff changeset
85 default_arg('type','char');
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
86 switch type
291
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
87 case{'c','char'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
88 [closure,penalty] = boundary_condition_char(obj,boundary);
293
2d604d16842c Works with varying coefficients and char boundary condition
Ylva Rydin <ylva.rydin@telia.com>
parents: 292
diff changeset
89 case{'general'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
90 [closure,penalty] = boundary_condition_general(obj,boundary,L);
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
91 otherwise
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
92 error('No such boundary condition')
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
93 end
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
94 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
95
946
706d1c2b4199 Raname opts to type in a bunch of interface methods
Jonatan Werpers <jonatan@werpers.com>
parents: 910
diff changeset
96 function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type)
706d1c2b4199 Raname opts to type in a bunch of interface methods
Jonatan Werpers <jonatan@werpers.com>
parents: 910
diff changeset
97 error('Not implemented');
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
98 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
99
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
100 function N = size(obj)
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
101 N = obj.m;
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
102 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
103
295
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
104 function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y)
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
105 params = obj.params;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
106
291
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
107 if isa(mat,'function_handle')
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
108 [rows,cols] = size(mat(params,0,0));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
109 matVec = mat(params,X',Y');
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
110 matVec = sparse(matVec);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
111 side = max(length(X),length(Y));
294
8ff6ec6249e8 "General" boundary conditions implemented
Ylva Rydin <ylva.rydin@telia.com>
parents: 293
diff changeset
112 else
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
113 matVec = mat;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
114 [rows,cols] = size(matVec);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
115 side = max(length(X),length(Y));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
116 cols = cols/side;
294
8ff6ec6249e8 "General" boundary conditions implemented
Ylva Rydin <ylva.rydin@telia.com>
parents: 293
diff changeset
117 end
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
118 ret = cell(rows,cols);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
119
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
120 for ii = 1:rows
294
8ff6ec6249e8 "General" boundary conditions implemented
Ylva Rydin <ylva.rydin@telia.com>
parents: 293
diff changeset
121 for jj=1:cols
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
122 ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side));
294
8ff6ec6249e8 "General" boundary conditions implemented
Ylva Rydin <ylva.rydin@telia.com>
parents: 293
diff changeset
123 end
8ff6ec6249e8 "General" boundary conditions implemented
Ylva Rydin <ylva.rydin@telia.com>
parents: 293
diff changeset
124 end
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
125 ret = cell2mat(ret);
294
8ff6ec6249e8 "General" boundary conditions implemented
Ylva Rydin <ylva.rydin@telia.com>
parents: 293
diff changeset
126 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
127
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
128 %Characteristic boundary conditions
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
129 function [closure, penalty] = boundary_condition_char(obj,boundary)
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
130 params = obj.params;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
131 x = obj.x;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
132 y = obj.y;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
133
291
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
134 switch boundary
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
135 case {'w','W','west'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
136 e_ = obj.e_w;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
137 mat = obj.A;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
138 boundPos = 'l';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
139 Hi = obj.Hxi;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
140 [V,Vi,D,signVec] = obj.matrixDiag(mat,x(1),y);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
141 side = max(length(y));
291
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
142 case {'e','E','east'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
143 e_ = obj.e_e;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
144 mat = obj.A;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
145 boundPos = 'r';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
146 Hi = obj.Hxi;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
147 [V,Vi,D,signVec] = obj.matrixDiag(mat,x(end),y);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
148 side = max(length(y));
291
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
149 case {'s','S','south'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
150 e_ = obj.e_s;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
151 mat = obj.B;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
152 boundPos = 'l';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
153 Hi = obj.Hyi;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
154 [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(1));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
155 side = max(length(x));
291
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
156 case {'n','N','north'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
157 e_ = obj.e_n;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
158 mat = obj.B;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
159 boundPos = 'r';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
160 Hi = obj.Hyi;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
161 [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(end));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
162 side = max(length(x));
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
163 end
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
164 pos = signVec(1);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
165 zeroval = signVec(2);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
166 neg = signVec(3);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
167
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
168 switch boundPos
294
8ff6ec6249e8 "General" boundary conditions implemented
Ylva Rydin <ylva.rydin@telia.com>
parents: 293
diff changeset
169 case {'l'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
170 tau = sparse(obj.n*side,pos);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
171 Vi_plus = Vi(1:pos,:);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
172 tau(1:pos,:) = -abs(D(1:pos,1:pos));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
173 closure = Hi*e_*V*tau*Vi_plus*e_';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
174 penalty = -Hi*e_*V*tau*Vi_plus;
294
8ff6ec6249e8 "General" boundary conditions implemented
Ylva Rydin <ylva.rydin@telia.com>
parents: 293
diff changeset
175 case {'r'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
176 tau = sparse(obj.n*side,neg);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
177 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
178 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
179 closure = Hi*e_*V*tau*Vi_minus*e_';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
180 penalty = -Hi*e_*V*tau*Vi_minus;
294
8ff6ec6249e8 "General" boundary conditions implemented
Ylva Rydin <ylva.rydin@telia.com>
parents: 293
diff changeset
181 end
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
182 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
183
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
184 % General boundary condition in the form Lu=g(x)
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
185 function [closure,penalty] = boundary_condition_general(obj,boundary,L)
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
186 params = obj.params;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
187 x = obj.x;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
188 y = obj.y;
997
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
189 e_ = obj.getBoundaryOperator('e', boundary);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
190
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
191 switch boundary
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
192 case {'w','W','west'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
193 mat = obj.A;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
194 boundPos = 'l';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
195 Hi = obj.Hxi;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
196 [V,Vi,D,signVec] = obj.matrixDiag(mat,x(1),y);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
197 L = obj.evaluateCoefficientMatrix(L,x(1),y);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
198 side = max(length(y));
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
199 case {'e','E','east'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
200 mat = obj.A;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
201 boundPos = 'r';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
202 Hi = obj.Hxi;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
203 [V,Vi,D,signVec] = obj.matrixDiag(mat,x(end),y);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
204 L = obj.evaluateCoefficientMatrix(L,x(end),y);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
205 side = max(length(y));
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
206 case {'s','S','south'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
207 mat = obj.B;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
208 boundPos = 'l';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
209 Hi = obj.Hyi;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
210 [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(1));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
211 L = obj.evaluateCoefficientMatrix(L,x,y(1));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
212 side = max(length(x));
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
213 case {'n','N','north'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
214 mat = obj.B;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
215 boundPos = 'r';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
216 Hi = obj.Hyi;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
217 [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(end));
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
218 L = obj.evaluateCoefficientMatrix(L,x,y(end));
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
219 side = max(length(x));
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
220 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
221
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
222 pos = signVec(1);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
223 zeroval = signVec(2);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
224 neg = signVec(3);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
225
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
226 switch boundPos
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
227 case {'l'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
228 tau = sparse(obj.n*side,pos);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
229 Vi_plus = Vi(1:pos,:);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
230 Vi_minus = Vi(pos+zeroval+1:obj.n*side,:);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
231 V_plus = V(:,1:pos);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
232 V_minus = V(:,(pos+zeroval)+1:obj.n*side);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
233
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
234 tau(1:pos,:) = -abs(D(1:pos,1:pos));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
235 R = -inv(L*V_plus)*(L*V_minus);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
236 closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
237 penalty = -Hi*e_*V*tau*inv(L*V_plus)*L;
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
238 case {'r'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
239 tau = sparse(obj.n*side,neg);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
240 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
241 Vi_plus = Vi(1:pos,:);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
242 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
243
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
244 V_plus = V(:,1:pos);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
245 V_minus = V(:,(pos+zeroval)+1:obj.n*side);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
246 R = -inv(L*V_minus)*(L*V_plus);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
247 closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
248 penalty = -Hi*e_*V*tau*inv(L*V_minus)*L;
294
8ff6ec6249e8 "General" boundary conditions implemented
Ylva Rydin <ylva.rydin@telia.com>
parents: 293
diff changeset
249 end
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
250 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
251
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
252 % Function that diagonalizes a symbolic matrix A as A=V*D*Vi
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
253 % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
254 % [d+ ]
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
255 % D = [ d0 ]
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
256 % [ d-]
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
257 % signVec is a vector specifying the number of possitive, zero and negative eigenvalues of D
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
258 function [V,Vi, D,signVec] = matrixDiag(obj,mat,x,y)
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
259 params = obj.params;
301
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 297
diff changeset
260 syms xs ys
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
261 [V, D]= eig(mat(params,xs,ys));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
262 Vi = inv(V);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
263 xs = x;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
264 ys = y;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
265
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
266 side = max(length(x),length(y));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
267 Dret = zeros(obj.n,side*obj.n);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
268 Vret = zeros(obj.n,side*obj.n);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
269 Viret = zeros(obj.n,side*obj.n);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
270
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
271 for ii = 1:obj.n
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
272 for jj = 1:obj.n
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
273 Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
274 Vret(jj,(ii-1)*side+1:side*ii) = eval(V(jj,ii));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
275 Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii));
291
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
276 end
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
277 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
278
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
279 D = sparse(Dret);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
280 V = sparse(Vret);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
281 Vi = sparse(Viret);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
282 V = obj.evaluateCoefficientMatrix(V,x,y);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
283 Vi = obj.evaluateCoefficientMatrix(Vi,x,y);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
284 D = obj.evaluateCoefficientMatrix(D,x,y);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
285 DD = diag(D);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
286
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
287 poseig = (DD>0);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
288 zeroeig = (DD==0);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
289 negeig = (DD<0);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
290
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
291 D = diag([DD(poseig); DD(zeroeig); DD(negeig)]);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
292 V = [V(:,poseig) V(:,zeroeig) V(:,negeig)];
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
293 Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
294 signVec = [sum(poseig),sum(zeroeig),sum(negeig)];
291
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
295 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
296
997
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
297 % Returns the boundary operator op for the boundary specified by the string boundary.
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
298 % op -- string or a cell array of strings
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
299 % boundary -- string
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
300 function varargout = getBoundaryOperator(obj, op, boundary)
1042
8d73fcdb07a5 Add asserts to boundary identifier inputs
Jonatan Werpers <jonatan@werpers.com>
parents: 997
diff changeset
301 assertIsMember(boundary, {'w', 'e', 's', 'n'})
997
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
302
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
303 if ~iscell(op)
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
304 op = {op};
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
305 end
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
306
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
307 for i = 1:numel(op)
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
308 switch op{i}
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
309 case 'e'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
310 switch boundary
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
311 case 'w'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
312 e = obj.e_w;
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
313 case 'e'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
314 e = obj.e_e;
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
315 case 's'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
316 e = obj.e_s;
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
317 case 'n'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
318 e = obj.e_n;
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
319 end
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
320 varargout{i} = e;
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
321 end
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
322 end
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
323 end
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
324
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
325 % Returns square boundary quadrature matrix, of dimension
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
326 % corresponding to the number of boundary points
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
327 %
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
328 % boundary -- string
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
329 function H_b = getBoundaryQuadrature(obj, boundary)
1042
8d73fcdb07a5 Add asserts to boundary identifier inputs
Jonatan Werpers <jonatan@werpers.com>
parents: 997
diff changeset
330 assertIsMember(boundary, {'w', 'e', 's', 'n'})
997
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
331
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
332 e = obj.getBoundaryOperator('e', boundary);
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
333
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
334 switch boundary
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
335 case 'w'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
336 H_b = inv(e'*obj.Hyi*e);
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
337 case 'e'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
338 H_b = inv(e'*obj.Hyi*e);
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
339 case 's'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
340 H_b = inv(e'*obj.Hxi*e);
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
341 case 'n'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
342 H_b = inv(e'*obj.Hxi*e);
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
343 end
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
344 end
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
345
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
346 end
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
347 end