Mercurial > repos > public > sbplib

annotate +scheme/Hypsyst3dCurve.m @ 1251:6424745e1b58 feature/volcano

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge in latest changes from default
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 20 Nov 2019 14:26:57 -0800
parents 0652b34f9f27
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
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
1226
0652b34f9f27 Add missing properties and methods to Hypsyst3D and Hypsyst3Dcurve. Note: These functions are NOT implemented properly and mearly work as place holders.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 946
diff changeset
6 grid %TODO: Abstract class requires a grid. Pass grid instead?
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
7 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
8 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
9
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
10 xi,eta,zeta
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
11 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
12
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 368
diff changeset
13 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
14 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
15
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
16 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
17
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
18 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
19 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
20 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
21 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
22
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 368
diff changeset
23 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
24
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
25 H % Discrete norm
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
26 % 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
27 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
28 Hxi,Heta,Hzeta
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
29 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
30 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
31 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
32 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
33 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
34
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
35
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
36 methods
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
37 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
38 xilim ={0 1};
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
39 etalim = {0 1};
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
40 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
41
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
42 if length(m) == 1
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
43 m = [m m m];
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
44 end
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
45 m_xi = m(1);
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
46 m_eta = m(2);
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
47 m_zeta = m(3);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
48 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
49 obj.params = params;
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
50 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
51
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
52 obj.m = m;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
53 obj.order = order;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
54 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
55
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
56 switch operator
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
57 case 'upwind'
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
58 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
59 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
60 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
61 case 'standard'
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
62 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
63 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
64 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
65 otherwise
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
66 error('Operator not available')
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
67 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
68
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
69 obj.xi = ops_xi.x;
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
70 obj.eta = ops_eta.x;
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
71 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
72
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 368
diff changeset
73 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
74 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
75 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
76
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
77
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
78 [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
79 obj.X = X;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
80 obj.Y = Y;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
81 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
82
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
83 I_n = eye(obj.n);
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
84 I_xi = speye(m_xi);
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
85 obj.I_xi = I_xi;
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
86 I_eta = speye(m_eta);
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
87 obj.I_eta = I_eta;
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
88 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
89 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
90
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
91 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
92
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
93 O_xi = ones(m_xi,1);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
94 O_eta = ones(m_eta,1);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
95 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
96
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
97
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
98 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
99 obj.Heta = ops_eta.H;
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
100 obj.Hzeta = ops_zeta.H;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
101 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
102
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
103 switch operator
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
104 case 'upwind'
366
7ada2db63268 Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 355
diff changeset
105 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
106 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
107 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
108 otherwise
366
7ada2db63268 Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 355
diff changeset
109 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
110 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
111 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
112 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
113
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
114 obj.A = A;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
115 obj.B = B;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
116 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
117
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
118 obj.X_xi = D1_xi*X;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
119 obj.X_eta = D1_eta*X;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
120 obj.X_zeta = D1_zeta*X;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
121 obj.Y_xi = D1_xi*Y;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
122 obj.Y_eta = D1_eta*Y;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
123 obj.Y_zeta = D1_zeta*Y;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
124 obj.Z_xi = D1_xi*Z;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
125 obj.Z_eta = D1_eta*Z;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
126 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
127
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
128 obj.Ahat = @transform_coefficient_matrix;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
129 obj.Bhat = @transform_coefficient_matrix;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
130 obj.Chat = @transform_coefficient_matrix;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
131 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
132
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
133 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
134 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
135 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
136
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
137 switch operator
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
138 case 'upwind'
366
7ada2db63268 Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 355
diff changeset
139 clear D1_xi D1_eta D1_zeta
368
53abf04f5e4e Made some small changes.
Ylva Rydin <ylva.rydin@telia.com>
parents: 366
diff changeset
140 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)))));
53abf04f5e4e Made some small changes.
Ylva Rydin <ylva.rydin@telia.com>
parents: 366
diff changeset
141 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)))));
53abf04f5e4e Made some small changes.
Ylva Rydin <ylva.rydin@telia.com>
parents: 366
diff changeset
142 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
143
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
144 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
145 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
146 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
147 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
148
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
149 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
150
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 368
diff changeset
151 obj.Aevaluated = [];
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
152 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
153 temp = Am*Dpxi;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 368
diff changeset
154 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
155 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
156
355
69b078cf8072 Upwind doed not work in the non curve-linear case.
Ylva Rydin <ylva.rydin@telia.com>
parents: 354
diff changeset
157 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
158 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
159 temp = Bp*Dmeta;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 368
diff changeset
160 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
161 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
162
355
69b078cf8072 Upwind doed not work in the non curve-linear case.
Ylva Rydin <ylva.rydin@telia.com>
parents: 354
diff changeset
163 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
164 obj.Bevaluated = [];
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
165 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
166 temp = Bm*Dpeta;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 368
diff changeset
167 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
168 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
169
355
69b078cf8072 Upwind doed not work in the non curve-linear case.
Ylva Rydin <ylva.rydin@telia.com>
parents: 354
diff changeset
170 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
171 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
172 temp = Cp*Dmzeta;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 368
diff changeset
173 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
174 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
175
366
7ada2db63268 Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 355
diff changeset
176 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
177 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
178 obj.Cevaluated = [];
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
179 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
180 temp = Cm*Dpzeta;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 368
diff changeset
181 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
182 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
183
366
7ada2db63268 Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 355
diff changeset
184 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
185 +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
186 +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
187 -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
188 -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
189 -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
190
366
7ada2db63268 Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 355
diff changeset
191 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
192 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
193
355
69b078cf8072 Upwind doed not work in the non curve-linear case.
Ylva Rydin <ylva.rydin@telia.com>
parents: 354
diff changeset
194 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
195
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 368
diff changeset
196 case 'standard'
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 368
diff changeset
197 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
198 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
199 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
200
366
7ada2db63268 Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 355
diff changeset
201 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
202 +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
203 +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
204 -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
205 -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
206 -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
207
366
7ada2db63268 Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 355
diff changeset
208 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
209 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
210
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
211 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
212 otherwise
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 368
diff changeset
213 error('Operator not supported')
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
214 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
215
366
7ada2db63268 Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 355
diff changeset
216 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
217 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
218 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
219
366
7ada2db63268 Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 355
diff changeset
220 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
221 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
222 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
223 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
224 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
225 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
226
366
7ada2db63268 Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 355
diff changeset
227 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
228 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
229 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
230 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
231 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
232 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
233
366
7ada2db63268 Tried to speed up the matrix set-up a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 355
diff changeset
234 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
235 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
236 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
237 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
238 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
239 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
240 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
241
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
242 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
243 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
244 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
245 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
246 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
247
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
248
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
249 % 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
250 % 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
251 % 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
252 % 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
253 % 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
254 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
255 default_arg('type','char');
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
256 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
257
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
258 switch type
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
259 case{'c','char'}
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
260 [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
261 case{'general'}
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
262 [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
263 otherwise
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
264 error('No such boundary condition')
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
265 end
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
266 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
267
946
706d1c2b4199 Raname opts to type in a bunch of interface methods
Jonatan Werpers <jonatan@werpers.com>
parents: 910
diff changeset
268 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
269 error('Not implemented');
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
270 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
271
1226
0652b34f9f27 Add missing properties and methods to Hypsyst3D and Hypsyst3Dcurve. Note: These functions are NOT implemented properly and mearly work as place holders.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 946
diff changeset
272 % TODO: Implement! This function should potentially replace boundary_matrices.
0652b34f9f27 Add missing properties and methods to Hypsyst3D and Hypsyst3Dcurve. Note: These functions are NOT implemented properly and mearly work as place holders.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 946
diff changeset
273 % Returns the boundary operator op for the boundary specified by the string boundary.
0652b34f9f27 Add missing properties and methods to Hypsyst3D and Hypsyst3Dcurve. Note: These functions are NOT implemented properly and mearly work as place holders.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 946
diff changeset
274 % op -- string or a cell array of strings
0652b34f9f27 Add missing properties and methods to Hypsyst3D and Hypsyst3Dcurve. Note: These functions are NOT implemented properly and mearly work as place holders.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 946
diff changeset
275 % boundary -- string
0652b34f9f27 Add missing properties and methods to Hypsyst3D and Hypsyst3Dcurve. Note: These functions are NOT implemented properly and mearly work as place holders.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 946
diff changeset
276 function varargout = getBoundaryOperator(obj, op, boundary)
0652b34f9f27 Add missing properties and methods to Hypsyst3D and Hypsyst3Dcurve. Note: These functions are NOT implemented properly and mearly work as place holders.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 946
diff changeset
277 error('Not implemented');
0652b34f9f27 Add missing properties and methods to Hypsyst3D and Hypsyst3Dcurve. Note: These functions are NOT implemented properly and mearly work as place holders.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 946
diff changeset
278 end
0652b34f9f27 Add missing properties and methods to Hypsyst3D and Hypsyst3Dcurve. Note: These functions are NOT implemented properly and mearly work as place holders.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 946
diff changeset
279
0652b34f9f27 Add missing properties and methods to Hypsyst3D and Hypsyst3Dcurve. Note: These functions are NOT implemented properly and mearly work as place holders.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 946
diff changeset
280 % TODO: Implement!
0652b34f9f27 Add missing properties and methods to Hypsyst3D and Hypsyst3Dcurve. Note: These functions are NOT implemented properly and mearly work as place holders.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 946
diff changeset
281 % Returns square boundary quadrature matrix, of dimension
0652b34f9f27 Add missing properties and methods to Hypsyst3D and Hypsyst3Dcurve. Note: These functions are NOT implemented properly and mearly work as place holders.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 946
diff changeset
282 % corresponding to the number of boundary points
0652b34f9f27 Add missing properties and methods to Hypsyst3D and Hypsyst3Dcurve. Note: These functions are NOT implemented properly and mearly work as place holders.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 946
diff changeset
283 %
0652b34f9f27 Add missing properties and methods to Hypsyst3D and Hypsyst3Dcurve. Note: These functions are NOT implemented properly and mearly work as place holders.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 946
diff changeset
284 % boundary -- string
0652b34f9f27 Add missing properties and methods to Hypsyst3D and Hypsyst3Dcurve. Note: These functions are NOT implemented properly and mearly work as place holders.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 946
diff changeset
285 function H = getBoundaryQuadrature(obj, boundary)
0652b34f9f27 Add missing properties and methods to Hypsyst3D and Hypsyst3Dcurve. Note: These functions are NOT implemented properly and mearly work as place holders.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 946
diff changeset
286 error('Not implemented');
0652b34f9f27 Add missing properties and methods to Hypsyst3D and Hypsyst3Dcurve. Note: These functions are NOT implemented properly and mearly work as place holders.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 946
diff changeset
287 end
0652b34f9f27 Add missing properties and methods to Hypsyst3D and Hypsyst3Dcurve. Note: These functions are NOT implemented properly and mearly work as place holders.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 946
diff changeset
288
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
289 function N = size(obj)
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
290 N = obj.m;
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
291 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
292
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 368
diff changeset
293 % 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
294 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
295 params = obj.params;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
296 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
297 if isa(mat,'function_handle')
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
298 [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
299 x_1 = kr(obj.onesN,x_1);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
300 x_2 = kr(obj.onesN,x_2);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
301 y_1 = kr(obj.onesN,y_1);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
302 y_2 = kr(obj.onesN,y_2);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
303 z_1 = kr(obj.onesN,z_1);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
304 z_2 = kr(obj.onesN,z_2);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
305 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
306 matVec = sparse(matVec);
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
307 else
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
308 matVec = mat;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
309 [rows,cols] = size(matVec);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
310 side = max(length(X),length(Y));
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
311 cols = cols/side;
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
312 end
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 368
diff changeset
313 matVec(abs(matVec)<10^(-10)) = 0;
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
314 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
315
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 368
diff changeset
316 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
317 for jj = 1:cols
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
318 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
319 end
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
320 end
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
321 ret = cell2mat(ret);
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
322 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
323
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
324 function [BM] = boundary_matrices(obj,boundary)
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
325 params = obj.params;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
326 BM.boundary = boundary;
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
327 switch boundary
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
328 case {'w','W','west'}
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
329 BM.e_ = obj.e_w;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
330 mat = obj.Ahat;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
331 BM.boundpos = 'l';
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
332 BM.Hi = obj.Hxii;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
333 BM.index = obj.index_w;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
334 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
335 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
336 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
337 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
338 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
339 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
340 case {'e','E','east'}
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
341 BM.e_ = obj.e_e;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
342 mat = obj.Ahat;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
343 BM.boundpos = 'r';
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
344 BM.Hi = obj.Hxii;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
345 BM.index = obj.index_e;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
346 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
347 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
348 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
349 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
350 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
351 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
352 case {'s','S','south'}
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
353 BM.e_ = obj.e_s;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
354 mat = obj.Bhat;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
355 BM.boundpos = 'l';
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
356 BM.Hi = obj.Hetai;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
357 BM.index = obj.index_s;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
358 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
359 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
360 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
361 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
362 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
363 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
364 case {'n','N','north'}
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
365 BM.e_ = obj.e_n;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
366 mat = obj.Bhat;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
367 BM.boundpos = 'r';
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
368 BM.Hi = obj.Hetai;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
369 BM.index = obj.index_n;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
370 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
371 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
372 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
373 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
374 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
375 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
376 case{'b','B','Bottom'}
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
377 BM.e_ = obj.e_b;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
378 mat = obj.Chat;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
379 BM.boundpos = 'l';
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
380 BM.Hi = obj.Hzetai;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
381 BM.index = obj.index_b;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
382 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
383 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
384 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
385 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
386 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
387 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
388 case{'t','T','Top'}
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
389 BM.e_ = obj.e_t;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
390 mat = obj.Chat;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
391 BM.boundpos = 'r';
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
392 BM.Hi = obj.Hzetai;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
393 BM.index = obj.index_t;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
394 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
395 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
396 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
397 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
398 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
399 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
400 end
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
401 [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
402 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
403 BM.side = sum(BM.index);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
404 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
405 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
406
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 368
diff changeset
407 % Characteristic boundary condition
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 368
diff changeset
408 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
409 side = BM.side;
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
410 pos = BM.pos;
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
411 neg = BM.neg;
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
412 zeroval = BM.zeroval;
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
413 V = BM.V;
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
414 Vi = BM.Vi;
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
415 Hi = BM.Hi;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
416 D = BM.D;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
417 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
418
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
419 switch BM.boundpos
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
420 case {'l'}
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
421 tau = sparse(obj.n*side,pos);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
422 Vi_plus = Vi(1:pos,:);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
423 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
424 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
425 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
426 case {'r'}
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
427 tau = sparse(obj.n*side,neg);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
428 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
429 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
430 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
431 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
432 end
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
433 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
434
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 368
diff changeset
435 % 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
436 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
437 side = BM.side;
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
438 pos = BM.pos;
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
439 neg = BM.neg;
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
440 zeroval = BM.zeroval;
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
441 V = BM.V;
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
442 Vi = BM.Vi;
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
443 Hi = BM.Hi;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
444 D = BM.D;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
445 e_ = BM.e_;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
446 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
447
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
448 switch BM.boundary
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
449 case{'b','B','bottom'}
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
450 Ji_vec = diag(obj.Ji);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
451 Ji = diag(Ji_vec(index));
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
452 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
453 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
454 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
455
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
456 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
457 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
458
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
459 switch BM.boundpos
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
460 case {'l'}
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
461 tau = sparse(obj.n*side,pos);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
462 Vi_plus = Vi(1:pos,:);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
463 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
464 V_plus = V(:,1:pos);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
465 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
466
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
467 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
468 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
469 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
470 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
471 case {'r'}
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
472 tau = sparse(obj.n*side,neg);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
473 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
474 Vi_plus = Vi(1:pos,:);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
475 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
476
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
477 V_plus = V(:,1:pos);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
478 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
479 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
480 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
481 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
482 end
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
483 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
484
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 368
diff changeset
485 % 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
486 % 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
487 % [d+ ]
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 368
diff changeset
488 % D = [ d0 ]
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 368
diff changeset
489 % [ d-]
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 368
diff changeset
490 % 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
491 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
492 params = obj.params;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
493 eps = 10^(-10);
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
494 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
495 syms x_1s
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
496 else
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
497 x_1s = 0;
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
498 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
499
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
500 if(sum(abs(x_2))>eps)
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
501 syms x_2s;
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
502 else
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
503 x_2s = 0;
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
504 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
505
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(y_1))>eps)
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
508 syms y_1s
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 y_1s = 0;
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
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
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
513 if(sum(abs(y_2))>eps)
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
514 syms y_2s;
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
515 else
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
516 y_2s = 0;
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
517 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
518
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
519 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
520 syms z_1s
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
521 else
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
522 z_1s = 0;
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
523 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
524
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
525 if(sum(abs(z_2))>eps)
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
526 syms z_2s;
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
527 else
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
528 z_2s = 0;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
529 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
530
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
531 syms xs ys zs
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
532 [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
533 Vi = inv(V);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
534 xs = x;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
535 ys = y;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
536 zs = z;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
537 x_1s = x_1;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
538 x_2s = x_2;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
539 y_1s = y_1;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
540 y_2s = y_2;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
541 z_1s = z_1;
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
542 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
543
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
544 side = max(length(x),length(y));
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
545 Dret = zeros(obj.n,side*obj.n);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
546 Vret = zeros(obj.n,side*obj.n);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
547 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
548
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
549 for ii=1:obj.n
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
550 for jj=1:obj.n
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
551 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
552 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
553 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
554 end
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
555 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
556
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
557 D = sparse(Dret);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
558 V = sparse(Vret);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
559 Vi = sparse(Viret);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
560 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
561 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
562 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
563 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
564
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
565 poseig = (DD>0);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
566 zeroeig = (DD==0);
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
567 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
568
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
569 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
570 V = [V(:,poseig) V(:,zeroeig) V(:,negeig)];
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
571 Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
572 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
573 end
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
574 end
354
dbac99d2c318 Removed inv(Vi) to save time
Ylva Rydin <ylva.rydin@telia.com>
parents: 351
diff changeset
575 end