Mercurial > repos > public > sbplib
annotate +scheme/Hypsyst3dCurve.m @ 1195:a4c00628a39d feature/rv
Add higher order approximations to BDFDerivative
author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
---|---|
date | Wed, 07 Aug 2019 13:27:36 +0200 |
parents | 706d1c2b4199 |
children | 0652b34f9f27 |
rev | line source |
---|---|
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
1 classdef Hypsyst3dCurve < scheme.Scheme |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
2 properties |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
3 m % Number of points in each direction, possibly a vector |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
4 n %size of system |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
5 h % Grid spacing |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
6 X, Y, Z% Values of x and y for each grid point |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
7 Yx, Zx, Xy, Zy, Xz, Yz %Grid values for boundary surfaces |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
8 |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
9 xi,eta,zeta |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
10 Xi, Eta, Zeta |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
11 |
369
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
12 Eta_xi, Zeta_xi, Xi_eta, Zeta_eta, Xi_zeta, Eta_zeta % Metric terms |
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
13 X_xi, X_eta, X_zeta,Y_xi,Y_eta,Y_zeta,Z_xi,Z_eta,Z_zeta % Metric terms |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
14 |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
15 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
|
16 |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
17 D % non-stabalized scheme operator |
369
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
18 Aevaluated, Bevaluated, Cevaluated, Eevaluated % Numeric Coeffiecient matrices |
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
19 Ahat, Bhat, Chat % Symbolic Transformed Coefficient matrices |
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
20 A, B, C, E % Symbolic coeffiecient 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
|
21 |
369
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
22 J, Ji % JAcobian and inverse Jacobian |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
23 |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
24 H % Discrete norm |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
25 % Norms in the x, y and z directions |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
26 Hxii,Hetai,Hzetai, Hzi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
27 Hxi,Heta,Hzeta |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
28 I_xi,I_eta,I_zeta, I_N,onesN |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
29 e_w, e_e, e_s, e_n, e_b, e_t |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
30 index_w, index_e,index_s,index_n, index_b, index_t |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
31 params %parameters for the coeficient matrice |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
32 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
|
33 |
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
34 |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
35 methods |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
36 function obj = Hypsyst3dCurve(m, order, A, B,C, E, params,ti,operator) |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
37 xilim ={0 1}; |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
38 etalim = {0 1}; |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
39 zetalim = {0 1}; |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
40 |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
41 if length(m) == 1 |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
42 m = [m m m]; |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
43 end |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
44 m_xi = m(1); |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
45 m_eta = m(2); |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
46 m_zeta = m(3); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
47 m_tot = m_xi*m_eta*m_zeta; |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
48 obj.params = params; |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
49 obj.n = length(A(obj,0,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
|
50 |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
51 obj.m = m; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
52 obj.order = order; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
53 obj.onesN = ones(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
|
54 |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
55 switch operator |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
56 case 'upwind' |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
57 ops_xi = sbp.D1Upwind(m_xi,xilim,order); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
58 ops_eta = sbp.D1Upwind(m_eta,etalim,order); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
59 ops_zeta = sbp.D1Upwind(m_zeta,zetalim,order); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
60 case 'standard' |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
61 ops_xi = sbp.D2Standard(m_xi,xilim,order); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
62 ops_eta = sbp.D2Standard(m_eta,etalim,order); |
366
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
63 ops_zeta = sbp.D2Standard(m_zeta,zetalim,order); |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
64 otherwise |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
65 error('Operator not available') |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
66 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
|
67 |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
68 obj.xi = ops_xi.x; |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
69 obj.eta = ops_eta.x; |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
70 obj.zeta = ops_zeta.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
|
71 |
369
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
72 obj.Xi = kr(obj.xi,ones(m_eta,1),ones(m_zeta,1)); |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
73 obj.Eta = kr(ones(m_xi,1),obj.eta,ones(m_zeta,1)); |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
74 obj.Zeta = kr(ones(m_xi,1),ones(m_eta,1),obj.zeta); |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
75 |
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
76 |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
77 [X,Y,Z] = ti.map(obj.Xi,obj.Eta,obj.Zeta); |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
78 obj.X = X; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
79 obj.Y = Y; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
80 obj.Z = Z; |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
81 |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
82 I_n = eye(obj.n); |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
83 I_xi = speye(m_xi); |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
84 obj.I_xi = I_xi; |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
85 I_eta = speye(m_eta); |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
86 obj.I_eta = I_eta; |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
87 I_zeta = speye(m_zeta); |
366
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
88 obj.I_zeta = I_zeta; |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
89 |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
90 I_N=kr(I_n,I_xi,I_eta,I_zeta); |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
91 |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
92 O_xi = ones(m_xi,1); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
93 O_eta = ones(m_eta,1); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
94 O_zeta = ones(m_zeta,1); |
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 |
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
96 |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
97 obj.Hxi = ops_xi.H; |
366
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
98 obj.Heta = ops_eta.H; |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
99 obj.Hzeta = ops_zeta.H; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
100 obj.h = [ops_xi.h ops_eta.h ops_zeta.h]; |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
101 |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
102 switch operator |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
103 case 'upwind' |
366
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
104 D1_xi = kr((ops_xi.Dp+ops_xi.Dm)/2, I_eta,I_zeta); |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
105 D1_eta = kr(I_xi, (ops_eta.Dp+ops_eta.Dm)/2,I_zeta); |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
106 D1_zeta = kr(I_xi, I_eta,(ops_zeta.Dp+ops_zeta.Dm)/2); |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
107 otherwise |
366
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
108 D1_xi = kr(ops_xi.D1, I_eta,I_zeta); |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
109 D1_eta = kr(I_xi, ops_eta.D1,I_zeta); |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
110 D1_zeta = kr(I_xi, I_eta,ops_zeta.D1); |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
111 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
|
112 |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
113 obj.A = A; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
114 obj.B = B; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
115 obj.C = C; |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
116 |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
117 obj.X_xi = D1_xi*X; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
118 obj.X_eta = D1_eta*X; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
119 obj.X_zeta = D1_zeta*X; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
120 obj.Y_xi = D1_xi*Y; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
121 obj.Y_eta = D1_eta*Y; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
122 obj.Y_zeta = D1_zeta*Y; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
123 obj.Z_xi = D1_xi*Z; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
124 obj.Z_eta = D1_eta*Z; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
125 obj.Z_zeta = D1_zeta*Z; |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
126 |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
127 obj.Ahat = @transform_coefficient_matrix; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
128 obj.Bhat = @transform_coefficient_matrix; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
129 obj.Chat = @transform_coefficient_matrix; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
130 obj.E = @(obj,x,y,z,~,~,~,~,~,~)E(obj,x,y,z); |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
131 |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
132 obj.Aevaluated = obj.evaluateCoefficientMatrix(obj.Ahat,obj.X, obj.Y,obj.Z, obj.X_eta,obj.X_zeta,obj.Y_eta,obj.Y_zeta,obj.Z_eta,obj.Z_zeta); |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
133 obj.Bevaluated = obj.evaluateCoefficientMatrix(obj.Bhat,obj.X, obj.Y,obj.Z, obj.X_zeta,obj.X_xi,obj.Y_zeta,obj.Y_xi,obj.Z_zeta,obj.Z_xi); |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
134 obj.Cevaluated = obj.evaluateCoefficientMatrix(obj.Chat,obj.X,obj.Y,obj.Z, obj.X_xi,obj.X_eta,obj.Y_xi,obj.Y_eta,obj.Z_xi,obj.Z_eta); |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
135 |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
136 switch operator |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
137 case 'upwind' |
366
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
138 clear D1_xi D1_eta D1_zeta |
368 | 139 alphaA = max(abs(eig(obj.Ahat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_eta(end),obj.X_zeta(end),obj.Y_eta(end),obj.Y_zeta(end),obj.Z_eta(end),obj.Z_zeta(end))))); |
140 alphaB = max(abs(eig(obj.Bhat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_zeta(end),obj.X_xi(end),obj.Y_zeta(end),obj.Y_xi(end),obj.Z_zeta(end),obj.Z_xi(end))))); | |
141 alphaC = max(abs(eig(obj.Chat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_xi(end),obj.X_eta(end),obj.Y_xi(end),obj.Y_eta(end),obj.Z_xi(end),obj.Z_eta(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
|
142 |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
143 Ap = (obj.Aevaluated+alphaA*I_N)/2; |
355
69b078cf8072
Upwind doed not work in the non curve-linear case.
Ylva Rydin <ylva.rydin@telia.com>
parents:
354
diff
changeset
|
144 Dmxi = kr(I_n, ops_xi.Dm, I_eta,I_zeta); |
369
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
145 diffSum = -Ap*Dmxi; |
355
69b078cf8072
Upwind doed not work in the non curve-linear case.
Ylva Rydin <ylva.rydin@telia.com>
parents:
354
diff
changeset
|
146 clear Ap Dmxi |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
147 |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
148 Am = (obj.Aevaluated-alphaA*I_N)/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
|
149 |
369
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
150 obj.Aevaluated = []; |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
151 Dpxi = kr(I_n, ops_xi.Dp, I_eta,I_zeta); |
369
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
152 temp = Am*Dpxi; |
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
153 diffSum = diffSum-temp; |
355
69b078cf8072
Upwind doed not work in the non curve-linear case.
Ylva Rydin <ylva.rydin@telia.com>
parents:
354
diff
changeset
|
154 clear Am Dpxi |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
155 |
355
69b078cf8072
Upwind doed not work in the non curve-linear case.
Ylva Rydin <ylva.rydin@telia.com>
parents:
354
diff
changeset
|
156 Bp = (obj.Bevaluated+alphaB*I_N)/2; |
69b078cf8072
Upwind doed not work in the non curve-linear case.
Ylva Rydin <ylva.rydin@telia.com>
parents:
354
diff
changeset
|
157 Dmeta = kr(I_n, I_xi, ops_eta.Dm,I_zeta); |
369
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
158 temp = Bp*Dmeta; |
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
159 diffSum = diffSum-temp; |
355
69b078cf8072
Upwind doed not work in the non curve-linear case.
Ylva Rydin <ylva.rydin@telia.com>
parents:
354
diff
changeset
|
160 clear Bp Dmeta |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
161 |
355
69b078cf8072
Upwind doed not work in the non curve-linear case.
Ylva Rydin <ylva.rydin@telia.com>
parents:
354
diff
changeset
|
162 Bm = (obj.Bevaluated-alphaB*I_N)/2; |
369
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
163 obj.Bevaluated = []; |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
164 Dpeta = kr(I_n, I_xi, ops_eta.Dp,I_zeta); |
369
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
165 temp = Bm*Dpeta; |
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
166 diffSum = diffSum-temp; |
366
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
167 clear Bm Dpeta |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
168 |
355
69b078cf8072
Upwind doed not work in the non curve-linear case.
Ylva Rydin <ylva.rydin@telia.com>
parents:
354
diff
changeset
|
169 Cp = (obj.Cevaluated+alphaC*I_N)/2; |
69b078cf8072
Upwind doed not work in the non curve-linear case.
Ylva Rydin <ylva.rydin@telia.com>
parents:
354
diff
changeset
|
170 Dmzeta = kr(I_n, I_xi, I_eta,ops_zeta.Dm); |
369
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
171 temp = Cp*Dmzeta; |
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
172 diffSum = diffSum-temp; |
355
69b078cf8072
Upwind doed not work in the non curve-linear case.
Ylva Rydin <ylva.rydin@telia.com>
parents:
354
diff
changeset
|
173 clear Cp Dmzeta |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
174 |
366
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
175 Cm = (obj.Cevaluated-alphaC*I_N)/2; |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
176 clear I_N |
369
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
177 obj.Cevaluated = []; |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
178 Dpzeta = kr(I_n, I_xi, I_eta,ops_zeta.Dp); |
369
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
179 temp = Cm*Dpzeta; |
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
180 diffSum = diffSum-temp; |
366
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
181 clear Cm Dpzeta temp |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
182 |
366
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
183 obj.J = obj.X_xi.*obj.Y_eta.*obj.Z_zeta... |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
184 +obj.X_zeta.*obj.Y_xi.*obj.Z_eta... |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
185 +obj.X_eta.*obj.Y_zeta.*obj.Z_xi... |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
186 -obj.X_xi.*obj.Y_zeta.*obj.Z_eta... |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
187 -obj.X_eta.*obj.Y_xi.*obj.Z_zeta... |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
188 -obj.X_zeta.*obj.Y_eta.*obj.Z_xi; |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
189 |
366
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
190 obj.Ji = kr(I_n,spdiags(1./obj.J,0,m_tot,m_tot)); |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
191 obj.Eevaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,obj.Z,[],[],[],[],[],[]); |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
192 |
355
69b078cf8072
Upwind doed not work in the non curve-linear case.
Ylva Rydin <ylva.rydin@telia.com>
parents:
354
diff
changeset
|
193 obj.D = obj.Ji*diffSum-obj.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
|
194 |
369
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
195 case 'standard' |
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
196 D1_xi = kr(I_n,D1_xi); |
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
197 D1_eta = kr(I_n,D1_eta); |
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
198 D1_zeta = kr(I_n,D1_zeta); |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
199 |
366
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
200 obj.J = obj.X_xi.*obj.Y_eta.*obj.Z_zeta... |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
201 +obj.X_zeta.*obj.Y_xi.*obj.Z_eta... |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
202 +obj.X_eta.*obj.Y_zeta.*obj.Z_xi... |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
203 -obj.X_xi.*obj.Y_zeta.*obj.Z_eta... |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
204 -obj.X_eta.*obj.Y_xi.*obj.Z_zeta... |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
205 -obj.X_zeta.*obj.Y_eta.*obj.Z_xi; |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
206 |
366
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
207 obj.Ji = kr(I_n,spdiags(1./obj.J,0,m_tot,m_tot)); |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
208 obj.Eevaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,obj.Z,[],[],[],[],[],[]); |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
209 |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
210 obj.D = obj.Ji*(-obj.Aevaluated*D1_xi-obj.Bevaluated*D1_eta -obj.Cevaluated*D1_zeta)-obj.Eevaluated; |
369
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
211 otherwise |
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
212 error('Operator not supported') |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
213 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
|
214 |
366
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
215 obj.Hxii = kr(I_n, ops_xi.HI, I_eta,I_zeta); |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
216 obj.Hetai = kr(I_n, I_xi, ops_eta.HI,I_zeta); |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
217 obj.Hzetai = kr(I_n, I_xi,I_eta, ops_zeta.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
|
218 |
366
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
219 obj.index_w = (kr(ops_xi.e_l, O_eta,O_zeta)==1); |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
220 obj.index_e = (kr(ops_xi.e_r, O_eta,O_zeta)==1); |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
221 obj.index_s = (kr(O_xi, ops_eta.e_l,O_zeta)==1); |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
222 obj.index_n = (kr(O_xi, ops_eta.e_r,O_zeta)==1); |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
223 obj.index_b = (kr(O_xi, O_eta, ops_zeta.e_l)==1); |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
224 obj.index_t = (kr(O_xi, O_eta, ops_zeta.e_r)==1); |
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 |
366
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
226 obj.e_w = kr(I_n, ops_xi.e_l, I_eta,I_zeta); |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
227 obj.e_e = kr(I_n, ops_xi.e_r, I_eta,I_zeta); |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
228 obj.e_s = kr(I_n, I_xi, ops_eta.e_l,I_zeta); |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
229 obj.e_n = kr(I_n, I_xi, ops_eta.e_r,I_zeta); |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
230 obj.e_b = kr(I_n, I_xi, I_eta, ops_zeta.e_l); |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
231 obj.e_t = kr(I_n, I_xi, I_eta, ops_zeta.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
|
232 |
366
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
233 obj.Eta_xi = kr(obj.eta,ones(m_xi,1)); |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
234 obj.Zeta_xi = kr(ones(m_eta,1),obj.zeta); |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
235 obj.Xi_eta = kr(obj.xi,ones(m_zeta,1)); |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
236 obj.Zeta_eta = kr(ones(m_xi,1),obj.zeta); |
7ada2db63268
Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
355
diff
changeset
|
237 obj.Xi_zeta = kr(obj.xi,ones(m_eta,1)); |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
238 obj.Eta_zeta = kr(ones(m_zeta,1),obj.eta); |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
239 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
|
240 |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
241 function [ret] = transform_coefficient_matrix(obj,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2) |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
242 ret = obj.A(obj,x,y,z).*(y_1.*z_2-z_1.*y_2); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
243 ret = ret+obj.B(obj,x,y,z).*(x_2.*z_1-x_1.*z_2); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
244 ret = ret+obj.C(obj,x,y,z).*(x_1.*y_2-x_2.*y_1); |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
245 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
|
246 |
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
247 |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
248 % Closure functions return the opertors applied to the own doamin to close the boundary |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
249 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
250 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
251 % type is a string specifying the type of boundary condition if there are several. |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
252 % data is a function returning the data that should be applied at the boundary. |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
253 function [closure, penalty] = boundary_condition(obj,boundary,type,L) |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
254 default_arg('type','char'); |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
255 BM = boundary_matrices(obj,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
|
256 |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
257 switch type |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
258 case{'c','char'} |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
259 [closure,penalty] = boundary_condition_char(obj,BM); |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
260 case{'general'} |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
261 [closure,penalty] = boundary_condition_general(obj,BM,boundary,L); |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
262 otherwise |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
263 error('No such boundary condition') |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
264 end |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
265 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
|
266 |
946
706d1c2b4199
Raname opts to type in a bunch of interface methods
Jonatan Werpers <jonatan@werpers.com>
parents:
910
diff
changeset
|
267 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
|
268 error('Not implemented'); |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
269 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
|
270 |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
271 function N = size(obj) |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
272 N = obj.m; |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
273 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
|
274 |
369
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
275 % Evaluates the symbolic Coeffiecient matrix mat |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
276 function [ret] = evaluateCoefficientMatrix(obj,mat, X, Y, Z , x_1 , x_2 , y_1 , y_2 , z_1 , z_2) |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
277 params = obj.params; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
278 side = max(length(X),length(Y)); |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
279 if isa(mat,'function_handle') |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
280 [rows,cols] = size(mat(obj,0,0,0,0,0,0,0,0,0)); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
281 x_1 = kr(obj.onesN,x_1); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
282 x_2 = kr(obj.onesN,x_2); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
283 y_1 = kr(obj.onesN,y_1); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
284 y_2 = kr(obj.onesN,y_2); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
285 z_1 = kr(obj.onesN,z_1); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
286 z_2 = kr(obj.onesN,z_2); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
287 matVec = mat(obj,X',Y',Z',x_1',x_2',y_1',y_2',z_1',z_2'); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
288 matVec = sparse(matVec); |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
289 else |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
290 matVec = mat; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
291 [rows,cols] = size(matVec); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
292 side = max(length(X),length(Y)); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
293 cols = cols/side; |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
294 end |
369
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
295 matVec(abs(matVec)<10^(-10)) = 0; |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
296 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
|
297 |
369
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
298 for ii = 1:rows |
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
299 for jj = 1:cols |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
300 ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side)); |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
301 end |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
302 end |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
303 ret = cell2mat(ret); |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
304 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
|
305 |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
306 function [BM] = boundary_matrices(obj,boundary) |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
307 params = obj.params; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
308 BM.boundary = boundary; |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
309 switch boundary |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
310 case {'w','W','west'} |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
311 BM.e_ = obj.e_w; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
312 mat = obj.Ahat; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
313 BM.boundpos = 'l'; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
314 BM.Hi = obj.Hxii; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
315 BM.index = obj.index_w; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
316 BM.x_1 = obj.X_eta(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
317 BM.x_2 = obj.X_zeta(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
318 BM.y_1 = obj.Y_eta(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
319 BM.y_2 = obj.Y_zeta(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
320 BM.z_1 = obj.Z_eta(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
321 BM.z_2 = obj.Z_zeta(BM.index); |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
322 case {'e','E','east'} |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
323 BM.e_ = obj.e_e; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
324 mat = obj.Ahat; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
325 BM.boundpos = 'r'; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
326 BM.Hi = obj.Hxii; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
327 BM.index = obj.index_e; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
328 BM.x_1 = obj.X_eta(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
329 BM.x_2 = obj.X_zeta(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
330 BM.y_1 = obj.Y_eta(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
331 BM.y_2 = obj.Y_zeta(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
332 BM.z_1 = obj.Z_eta(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
333 BM.z_2 = obj.Z_zeta(BM.index); |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
334 case {'s','S','south'} |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
335 BM.e_ = obj.e_s; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
336 mat = obj.Bhat; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
337 BM.boundpos = 'l'; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
338 BM.Hi = obj.Hetai; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
339 BM.index = obj.index_s; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
340 BM.x_1 = obj.X_zeta(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
341 BM.x_2 = obj.X_xi(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
342 BM.y_1 = obj.Y_zeta(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
343 BM.y_2 = obj.Y_xi(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
344 BM.z_1 = obj.Z_zeta(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
345 BM.z_2 = obj.Z_xi(BM.index); |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
346 case {'n','N','north'} |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
347 BM.e_ = obj.e_n; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
348 mat = obj.Bhat; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
349 BM.boundpos = 'r'; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
350 BM.Hi = obj.Hetai; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
351 BM.index = obj.index_n; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
352 BM.x_1 = obj.X_zeta(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
353 BM.x_2 = obj.X_xi(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
354 BM.y_1 = obj.Y_zeta(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
355 BM.y_2 = obj.Y_xi(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
356 BM.z_1 = obj.Z_zeta(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
357 BM.z_2 = obj.Z_xi(BM.index); |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
358 case{'b','B','Bottom'} |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
359 BM.e_ = obj.e_b; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
360 mat = obj.Chat; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
361 BM.boundpos = 'l'; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
362 BM.Hi = obj.Hzetai; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
363 BM.index = obj.index_b; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
364 BM.x_1 = obj.X_xi(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
365 BM.x_2 = obj.X_eta(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
366 BM.y_1 = obj.Y_xi(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
367 BM.y_2 = obj.Y_eta(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
368 BM.z_1 = obj.Z_xi(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
369 BM.z_2 = obj.Z_eta(BM.index); |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
370 case{'t','T','Top'} |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
371 BM.e_ = obj.e_t; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
372 mat = obj.Chat; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
373 BM.boundpos = 'r'; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
374 BM.Hi = obj.Hzetai; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
375 BM.index = obj.index_t; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
376 BM.x_1 = obj.X_xi(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
377 BM.x_2 = obj.X_eta(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
378 BM.y_1 = obj.Y_xi(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
379 BM.y_2 = obj.Y_eta(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
380 BM.z_1 = obj.Z_xi(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
381 BM.z_2 = obj.Z_eta(BM.index); |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
382 end |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
383 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.X(BM.index),obj.Y(BM.index),obj.Z(BM.index),... |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
384 BM.x_1,BM.x_2,BM.y_1,BM.y_2,BM.z_1,BM.z_2); |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
385 BM.side = sum(BM.index); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
386 BM.pos = signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3); |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
387 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
|
388 |
369
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
389 % Characteristic boundary condition |
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
390 function [closure, penalty] = boundary_condition_char(obj,BM) |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
391 side = BM.side; |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
392 pos = BM.pos; |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
393 neg = BM.neg; |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
394 zeroval = BM.zeroval; |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
395 V = BM.V; |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
396 Vi = BM.Vi; |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
397 Hi = BM.Hi; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
398 D = BM.D; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
399 e_ = BM.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
|
400 |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
401 switch BM.boundpos |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
402 case {'l'} |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
403 tau = sparse(obj.n*side,pos); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
404 Vi_plus = Vi(1:pos,:); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
405 tau(1:pos,:) = -abs(D(1:pos,1:pos)); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
406 closure = Hi*e_*V*tau*Vi_plus*e_'; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
407 penalty = -Hi*e_*V*tau*Vi_plus; |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
408 case {'r'} |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
409 tau = sparse(obj.n*side,neg); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
410 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
411 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
412 closure = Hi*e_*V*tau*Vi_minus*e_'; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
413 penalty = -Hi*e_*V*tau*Vi_minus; |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
414 end |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
415 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
|
416 |
369
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
417 % 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:
368
diff
changeset
|
418 function [closure,penalty] = boundary_condition_general(obj,BM,boundary,L) |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
419 side = BM.side; |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
420 pos = BM.pos; |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
421 neg = BM.neg; |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
422 zeroval = BM.zeroval; |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
423 V = BM.V; |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
424 Vi = BM.Vi; |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
425 Hi = BM.Hi; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
426 D = BM.D; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
427 e_ = BM.e_; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
428 index = BM.index; |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
429 |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
430 switch BM.boundary |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
431 case{'b','B','bottom'} |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
432 Ji_vec = diag(obj.Ji); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
433 Ji = diag(Ji_vec(index)); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
434 Zeta_x = Ji*(obj.Y_xi(index).*obj.Z_eta(index)-obj.Z_xi(index).*obj.Y_eta(index)); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
435 Zeta_y = Ji*(obj.X_eta(index).*obj.Z_xi(index)-obj.X_xi(index).*obj.Z_eta(index)); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
436 Zeta_z = Ji*(obj.X_xi(index).*obj.Y_eta(index)-obj.Y_xi(index).*obj.X_eta(index)); |
905
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
437 |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
438 L = obj.evaluateCoefficientMatrix(L,Zeta_x,Zeta_y,Zeta_z,[],[],[],[],[],[]); |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
439 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
|
440 |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
441 switch BM.boundpos |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
442 case {'l'} |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
443 tau = sparse(obj.n*side,pos); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
444 Vi_plus = Vi(1:pos,:); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
445 Vi_minus = Vi(pos+zeroval+1:obj.n*side,:); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
446 V_plus = V(:,1:pos); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
447 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
|
448 |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
449 tau(1:pos,:) = -abs(D(1:pos,1:pos)); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
450 R = -inv(L*V_plus)*(L*V_minus); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
451 closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
452 penalty = -Hi*e_*V*tau*inv(L*V_plus)*L; |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
453 case {'r'} |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
454 tau = sparse(obj.n*side,neg); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
455 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
456 Vi_plus = Vi(1:pos,:); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
457 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
|
458 |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
459 V_plus = V(:,1:pos); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
460 V_minus = V(:,(pos+zeroval)+1:obj.n*side); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
461 R = -inv(L*V_minus)*(L*V_plus); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
462 closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_'; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
463 penalty = -Hi*e_*V*tau*inv(L*V_minus)*L; |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
464 end |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
465 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
|
466 |
369
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
467 % 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:
368
diff
changeset
|
468 % 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:
368
diff
changeset
|
469 % [d+ ] |
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
470 % D = [ d0 ] |
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
471 % [ d-] |
9d1fc984f40d
Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents:
368
diff
changeset
|
472 % signVec is a vector specifying the number of possitive, zero and negative eigenvalues of D |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
473 function [V,Vi, D,signVec] = matrixDiag(obj,mat,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2) |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
474 params = obj.params; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
475 eps = 10^(-10); |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
476 if(sum(abs(x_1))>eps) |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
477 syms x_1s |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
478 else |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
479 x_1s = 0; |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
480 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
|
481 |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
482 if(sum(abs(x_2))>eps) |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
483 syms x_2s; |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
484 else |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
485 x_2s = 0; |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
486 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
|
487 |
459eeb99130f
Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents:
369
diff
changeset
|
488 |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
489 if(sum(abs(y_1))>eps) |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
490 syms y_1s |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
491 else |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
492 y_1s = 0; |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
493 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
|
494 |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
495 if(sum(abs(y_2))>eps) |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
496 syms y_2s; |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
497 else |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
498 y_2s = 0; |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
499 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
|
500 |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
501 if(sum(abs(z_1))>eps) |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
502 syms z_1s |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
503 else |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
504 z_1s = 0; |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
505 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
|
506 |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
507 if(sum(abs(z_2))>eps) |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
508 syms z_2s; |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
509 else |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
510 z_2s = 0; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
511 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
|
512 |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
513 syms xs ys zs |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
514 [V, D] = eig(mat(obj,xs,ys,zs,x_1s,x_2s,y_1s,y_2s,z_1s,z_2s)); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
515 Vi = inv(V); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
516 xs = x; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
517 ys = y; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
518 zs = z; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
519 x_1s = x_1; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
520 x_2s = x_2; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
521 y_1s = y_1; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
522 y_2s = y_2; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
523 z_1s = z_1; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
524 z_2s = z_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
|
525 |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
526 side = max(length(x),length(y)); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
527 Dret = zeros(obj.n,side*obj.n); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
528 Vret = zeros(obj.n,side*obj.n); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
529 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
|
530 |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
531 for ii=1:obj.n |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
532 for jj=1:obj.n |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
533 Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii)); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
534 Vret(jj,(ii-1)*side+1:side*ii) = eval(V(jj,ii)); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
535 Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii)); |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
536 end |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
537 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
|
538 |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
539 D = sparse(Dret); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
540 V = sparse(Vret); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
541 Vi = sparse(Viret); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
542 V = obj.evaluateCoefficientMatrix(V,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
543 D = obj.evaluateCoefficientMatrix(D,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
544 Vi = obj.evaluateCoefficientMatrix(Vi,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
545 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
|
546 |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
547 poseig = (DD>0); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
548 zeroeig = (DD==0); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
549 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
|
550 |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
551 D = diag([DD(poseig); DD(zeroeig); DD(negeig)]); |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
552 V = [V(:,poseig) V(:,zeroeig) V(:,negeig)]; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
553 Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)]; |
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
554 signVec = [sum(poseig),sum(zeroeig),sum(negeig)]; |
350
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
555 end |
5d5652fe826a
A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff
changeset
|
556 end |
354
dbac99d2c318
Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents:
351
diff
changeset
|
557 end |