Mercurial > repos > public > sbplib

annotate +scheme/Hypsyst2dCurve.m @ 366:7ada2db63268 feature/hypsyst

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Tried to speed up the matrix set-up a little
author Ylva Rydin <ylva.rydin@telia.com>
date Tue, 13 Dec 2016 09:28:33 +0100
parents 9b3d7fc61a36
children 9d1fc984f40d
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
1 classdef Hypsyst2dCurve < scheme.Scheme
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
2 properties
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
3 m % Number of points in each direction, possibly a vector
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
4 n %size of system
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
5 h % Grid spacing
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
6 X,Y % Values of x and y for each grid point
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
7
301
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
8 J, Ji %Jacobaian and inverse Jacobian
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
9 xi,eta
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
10 Xi,Eta
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
11
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
12 A,B
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
13 X_eta, Y_eta
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
14 X_xi,Y_xi
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
15 order % Order accuracy for the approximation
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
16
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
17 D % non-stabalized scheme operator
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
18 Ahat, Bhat, E
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
19
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
20 H % Discrete norm
301
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
21 Hxii,Hetai % Kroneckerd norms in xi and eta.
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
22 I_xi,I_eta, I_N, onesN
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
23 e_w, e_e, e_s, e_n
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
24 index_w, index_e,index_s,index_n
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
25 params %parameters for the coeficient matrice
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
26 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
27
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
28
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
29 methods
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
30 function obj = Hypsyst2dCurve(m, order, A, B, E, params, ti)
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
31 default_arg('E', [])
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
32 xilim = {0 1};
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
33 etalim = {0 1};
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
34
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
35 if length(m) == 1
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
36 m = [m m];
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
37 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
38 obj.params = params;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
39 obj.A=A;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
40 obj.B=B;
301
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
41
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
42 obj.Ahat=@(params,x,y,x_eta,y_eta)(A(params,x,y).*y_eta-B(params,x,y).*x_eta);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
43 obj.Bhat=@(params,x,y,x_xi,y_xi)(B(params,x,y).*x_xi-A(params,x,y).*y_xi);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
44 obj.E=@(params,x,y,~,~)E(params,x,y);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
45
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
46 m_xi = m(1);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
47 m_eta = m(2);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
48 m_tot=m_xi*m_eta;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
49
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
50 ops_xi = sbp.D2Standard(m_xi,xilim,order);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
51 ops_eta = sbp.D2Standard(m_eta,etalim,order);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
52
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
53 obj.xi = ops_xi.x;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
54 obj.eta = ops_eta.x;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
55
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
56 obj.Xi = kr(obj.xi,ones(m_eta,1));
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
57 obj.Eta = kr(ones(m_xi,1),obj.eta);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
58
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
59 obj.n = length(A(obj.params,0,0));
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
60 obj.onesN=ones(obj.n);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
61
301
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
62 obj.index_w=1:m_eta;
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
63 obj.index_e=(m_tot-m_eta+1):m_tot;
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
64 obj.index_s=1:m_eta:(m_tot-m_eta+1);
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
65 obj.index_n=(m_eta):m_eta:m_tot;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
66
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
67 I_n = eye(obj.n);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
68 I_xi = speye(m_xi);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
69 obj.I_xi = I_xi;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
70 I_eta = speye(m_eta);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
71 obj.I_eta = I_eta;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
72
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
73 D1_xi = kr(I_n, ops_xi.D1, I_eta);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
74 obj.Hxii = kr(I_n, ops_xi.HI, I_eta);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
75 D1_eta = kr(I_n, I_xi, ops_eta.D1);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
76 obj.Hetai = kr(I_n, I_xi, ops_eta.HI);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
77
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
78 obj.e_w = kr(I_n, ops_xi.e_l, I_eta);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
79 obj.e_e = kr(I_n, ops_xi.e_r, I_eta);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
80 obj.e_s = kr(I_n, I_xi, ops_eta.e_l);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
81 obj.e_n = kr(I_n, I_xi, ops_eta.e_r);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
82
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
83 [X,Y] = ti.map(obj.xi,obj.eta);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
84
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
85 [x_xi,x_eta] = gridDerivatives(X,ops_xi.D1,ops_eta.D1);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
86 [y_xi,y_eta] = gridDerivatives(Y,ops_xi.D1, ops_eta.D1);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
87
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
88 obj.X=reshape(X,m_tot,1);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
89 obj.Y=reshape(Y,m_tot,1);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
90 obj.X_xi=reshape(x_xi,m_tot,1);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
91 obj.Y_xi=reshape(y_xi,m_tot,1);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
92 obj.X_eta=reshape(x_eta,m_tot,1);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
93 obj.Y_eta=reshape(y_eta,m_tot,1);
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
94
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
95 Ahat_evaluated = obj.evaluateCoefficientMatrix(obj.Ahat, obj.X, obj.Y,obj.X_eta,obj.Y_eta);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
96 Bhat_evaluated = obj.evaluateCoefficientMatrix(obj.Bhat, obj.X, obj.Y,obj.X_xi,obj.Y_xi);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
97 E_evaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,[],[]);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
98
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
99 obj.m=m;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
100 obj.h=[ops_xi.h ops_eta.h];
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
101 obj.order=order;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
102 obj.J=obj.X_xi.*obj.Y_eta - obj.X_eta.*obj.Y_xi;
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
103 obj.Ji =kr(I_n,spdiags(1./obj.J,0,m_tot,m_tot));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
104
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
105 obj.D=obj.Ji*(-Ahat_evaluated*D1_xi-Bhat_evaluated*D1_eta)-E_evaluated;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
106
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
107 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
108
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
109 % Closure functions return the opertors applied to the own doamin to close the boundary
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
110 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
111 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
112 % type is a string specifying the type of boundary condition if there are several.
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
113 % data is a function returning the data that should be applied at the boundary.
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
114 function [closure, penalty] = boundary_condition(obj,boundary,type,L)
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
115 default_arg('type','char');
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
116 switch type
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
117 case{'c','char'}
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
118 [closure,penalty]=boundary_condition_char(obj,boundary);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
119 case{'general'}
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
120 [closure,penalty]=boundary_condition_general(obj,boundary,L);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
121 otherwise
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
122 error('No such boundary condition')
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
123 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
124 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
125
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
126 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
127 error('An interface function does not exist yet');
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
128 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
129
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
130 function N = size(obj)
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
131 N = obj.m;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
132 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
133
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
134 function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y,x_,y_)
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
135 params=obj.params;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
136
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
137 if isa(mat,'function_handle')
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
138 [rows,cols]=size(mat(params,0,0,0,0));
301
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
139 x_=kr(obj.onesN,x_);
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
140 y_=kr(obj.onesN,y_);
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
141 matVec=mat(params,X',Y',x_',y_');
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
142 matVec=sparse(matVec);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
143 side=max(length(X),length(Y));
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
144 else
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
145 matVec=mat;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
146 [rows,cols]=size(matVec);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
147 side=max(length(X),length(Y));
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
148 cols=cols/side;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
149 end
352
9b3d7fc61a36 Changed operator in Utux
Ylva Rydin <ylva.rydin@telia.com>
parents: 350
diff changeset
150 ret=cell(rows,cols);
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
151
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
152 for ii=1:rows
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
153 for jj=1:cols
352
9b3d7fc61a36 Changed operator in Utux
Ylva Rydin <ylva.rydin@telia.com>
parents: 350
diff changeset
154 ret{ii,jj}=diag(matVec(ii,(jj-1)*side+1:jj*side));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
155 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
156 end
352
9b3d7fc61a36 Changed operator in Utux
Ylva Rydin <ylva.rydin@telia.com>
parents: 350
diff changeset
157 ret=cell2mat(ret);
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
158 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
159
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
160
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
161 function [closure, penalty]=boundary_condition_char(obj,boundary)
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
162 params=obj.params;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
163 X=obj.X; Y=obj.Y;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
164 xi=obj.xi; eta=obj.eta;
301
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
165
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
166
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
167 switch boundary
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
168 case {'w','W','west'}
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
169 e_=obj.e_w;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
170 mat=obj.Ahat;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
171 boundPos='l';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
172 Hi=obj.Hxii;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
173 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_w),Y(obj.index_w),obj.X_eta(obj.index_w),obj.Y_eta(obj.index_w));
301
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
174 side=max(length(eta));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
175 case {'e','E','east'}
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
176 e_=obj.e_e;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
177 mat=obj.Ahat;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
178 boundPos='r';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
179 Hi=obj.Hxii;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
180 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_e),Y(obj.index_e),obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e));
301
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
181 side=max(length(eta));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
182 case {'s','S','south'}
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
183 e_=obj.e_s;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
184 mat=obj.Bhat;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
185 boundPos='l';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
186 Hi=obj.Hetai;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
187 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_s),Y(obj.index_s),obj.X_xi(obj.index_s),obj.Y_xi(obj.index_s));
301
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
188 side=max(length(xi));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
189 case {'n','N','north'}
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
190 e_=obj.e_n;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
191 mat=obj.Bhat;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
192 boundPos='r';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
193 Hi=obj.Hetai;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
194 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_n),Y(obj.index_n),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n));
301
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
195 side=max(length(xi));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
196 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
197
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
198 pos=signVec(1); zeroval=signVec(2); neg=signVec(3);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
199
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
200 switch boundPos
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
201 case {'l'}
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
202 tau=sparse(obj.n*side,pos);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
203 Vi_plus=Vi(1:pos,:);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
204 tau(1:pos,:)=-abs(D(1:pos,1:pos));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
205 closure=Hi*e_*V*tau*Vi_plus*e_';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
206 penalty=-Hi*e_*V*tau*Vi_plus;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
207 case {'r'}
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
208 tau=sparse(obj.n*side,neg);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
209 tau((pos+zeroval)+1:obj.n*side,:)=-abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
210 Vi_minus=Vi((pos+zeroval)+1:obj.n*side,:);
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
211 closure=Hi*e_*V*tau*Vi_minus*e_';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
212 penalty=-Hi*e_*V*tau*Vi_minus;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
213 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
214 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
215
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
216
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
217 function [closure,penalty]=boundary_condition_general(obj,boundary,L)
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
218 params=obj.params;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
219 X=obj.X; Y=obj.Y;
301
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
220 xi=obj.xi; eta=obj.eta;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
221
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
222 switch boundary
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
223 case {'w','W','west'}
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
224 e_=obj.e_w;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
225 mat=obj.Ahat;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
226 boundPos='l';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
227 Hi=obj.Hxii;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
228 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_w),Y(obj.index_w),obj.X_eta(obj.index_w),obj.Y_eta(obj.index_w));
301
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
229
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
230 Ji_vec=diag(obj.Ji);
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
231 Ji=diag(Ji_vec(obj.index_w));
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
232 xi_x=Ji*obj.Y_eta(obj.index_w);
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
233 xi_y=-Ji*obj.X_eta(obj.index_w);
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
234 L=obj.evaluateCoefficientMatrix(L,xi_x,xi_y,[],[]);
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
235 side=max(length(eta));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
236 case {'e','E','east'}
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
237 e_=obj.e_e;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
238 mat=obj.Ahat;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
239 boundPos='r';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
240 Hi=obj.Hxii;
301
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
241 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_e),Y(obj.index_e),obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e));
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
242
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
243 Ji_vec=diag(obj.Ji);
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
244 Ji=diag(Ji_vec(obj.index_e));
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
245 xi_x=Ji*obj.Y_eta(obj.index_e);
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
246 xi_y=-Ji*obj.X_eta(obj.index_e);
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
247 L=obj.evaluateCoefficientMatrix(L,-xi_x,-xi_y,[],[]);
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
248 side=max(length(eta));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
249 case {'s','S','south'}
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
250 e_=obj.e_s;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
251 mat=obj.Bhat;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
252 boundPos='l';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
253 Hi=obj.Hetai;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
254 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_s),Y(obj.index_s),obj.X_xi(obj.index_s),obj.Y_xi(obj.index_s));
301
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
255
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
256 Ji_vec=diag(obj.Ji);
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
257 Ji=diag(Ji_vec(obj.index_s));
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
258 eta_x=Ji*obj.Y_xi(obj.index_s);
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
259 eta_y=-Ji*obj.X_xi(obj.index_s);
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
260 L=obj.evaluateCoefficientMatrix(L,eta_x,eta_y,[],[]);
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
261 side=max(length(xi));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
262 case {'n','N','north'}
352
9b3d7fc61a36 Changed operator in Utux
Ylva Rydin <ylva.rydin@telia.com>
parents: 350
diff changeset
263 e_=obj.e_n;
350
5d5652fe826a A commit before I try resolving the performance issues
Ylva Rydin <ylva.rydin@telia.com>
parents: 301
diff changeset
264
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
265 mat=obj.Bhat;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
266 boundPos='r';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
267 Hi=obj.Hetai;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
268 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_n),Y(obj.index_n),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n));
301
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
269
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
270 Ji_vec=diag(obj.Ji);
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
271 Ji=diag(Ji_vec(obj.index_n));
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
272 eta_x=Ji*obj.Y_xi(obj.index_n);
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
273 eta_y=-Ji*obj.X_xi(obj.index_n);
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
274 L=obj.evaluateCoefficientMatrix(L,-eta_x,-eta_y,[],[]);
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
275
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
276 side=max(length(xi));
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
277
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
278 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
279
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
280 pos=signVec(1); zeroval=signVec(2); neg=signVec(3);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
281
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
282 switch boundPos
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
283 case {'l'}
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
284 tau=sparse(obj.n*side,pos);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
285 Vi_plus=Vi(1:pos,:);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
286 Vi_minus=Vi(pos+1:obj.n*side,:);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
287 V_plus=V(:,1:pos);
301
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
288 V_minus=V(:,(pos)+1:obj.n*side);
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
289
301
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
290 tau(1:pos,:)=-abs(D(1:pos,1:pos));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
291 R=-inv(L*V_plus)*(L*V_minus);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
292 closure=Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
293 penalty=-Hi*e_*V*tau*inv(L*V_plus)*L;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
294 case {'r'}
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
295 tau=sparse(obj.n*side,neg);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
296 tau((pos+zeroval)+1:obj.n*side,:)=-abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
297 Vi_plus=Vi(1:pos,:);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
298 Vi_minus=Vi((pos+zeroval)+1:obj.n*side,:);
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
299
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
300 V_plus=V(:,1:pos);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
301 V_minus=V(:,(pos+zeroval)+1:obj.n*side);
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
302 R=-inv(L*V_minus)*(L*V_plus);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
303 closure=Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
304 penalty=-Hi*e_*V*tau*inv(L*V_minus)*L;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
305 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
306 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
307
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
308 function [V,Vi, D,signVec]=matrixDiag(obj,mat,x,y,x_,y_)
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
309 params=obj.params;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
310 syms xs ys
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
311 if(sum(abs(x_))~=0)
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
312 syms xs_
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
313 else
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
314 xs_=0;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
315 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
316
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
317 if(sum(abs(y_))~=0)
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
318 syms ys_;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
319 else
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
320 ys_=0;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
321 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
322
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
323 [V, D]=eig(mat(params,xs,ys,xs_,ys_));
352
9b3d7fc61a36 Changed operator in Utux
Ylva Rydin <ylva.rydin@telia.com>
parents: 350
diff changeset
324 Vi=inv(V);
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
325 syms xs ys xs_ ys_
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
326
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
327 xs=x;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
328 ys=y;
299
4d8d6eb0c116 A commit before I change the boundary treatment
Ylva Rydin <ylva.rydin@telia.com>
parents: 298
diff changeset
329 xs_=x_;
4d8d6eb0c116 A commit before I change the boundary treatment
Ylva Rydin <ylva.rydin@telia.com>
parents: 298
diff changeset
330 ys_=y_;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
331
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
332 side=max(length(x),length(y));
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
333 Dret=zeros(obj.n,side*obj.n);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
334 Vret=zeros(obj.n,side*obj.n);
352
9b3d7fc61a36 Changed operator in Utux
Ylva Rydin <ylva.rydin@telia.com>
parents: 350
diff changeset
335 Viret=zeros(obj.n,side*obj.n);
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
336 for ii=1:obj.n
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
337 for jj=1:obj.n
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
338 Dret(jj,(ii-1)*side+1:side*ii)=eval(D(jj,ii));
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
339 Vret(jj,(ii-1)*side+1:side*ii)=eval(V(jj,ii));
352
9b3d7fc61a36 Changed operator in Utux
Ylva Rydin <ylva.rydin@telia.com>
parents: 350
diff changeset
340 Viret(jj,(ii-1)*side+1:side*ii)=eval(Vi(jj,ii));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
341 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
342 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
343
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
344 D=sparse(Dret);
352
9b3d7fc61a36 Changed operator in Utux
Ylva Rydin <ylva.rydin@telia.com>
parents: 350
diff changeset
345 V=sparse(Vret);
9b3d7fc61a36 Changed operator in Utux
Ylva Rydin <ylva.rydin@telia.com>
parents: 350
diff changeset
346 Vi=sparse(Viret);
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
347 V=obj.evaluateCoefficientMatrix(V,x,y,x_,y_);
352
9b3d7fc61a36 Changed operator in Utux
Ylva Rydin <ylva.rydin@telia.com>
parents: 350
diff changeset
348 D=obj.evaluateCoefficientMatrix(D,x,y,x_,y_);
9b3d7fc61a36 Changed operator in Utux
Ylva Rydin <ylva.rydin@telia.com>
parents: 350
diff changeset
349 Vi=obj.evaluateCoefficientMatrix(Vi,x,y,x_,y_);
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
350 DD=diag(D);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
351
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
352 poseig=(DD>0);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
353 zeroeig=(DD==0);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
354 negeig=(DD<0);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
355
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
356 D=diag([DD(poseig); DD(zeroeig); DD(negeig)]);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
357 V=[V(:,poseig) V(:,zeroeig) V(:,negeig)];
352
9b3d7fc61a36 Changed operator in Utux
Ylva Rydin <ylva.rydin@telia.com>
parents: 350
diff changeset
358 Vi=[Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
359 signVec=[sum(poseig),sum(zeroeig),sum(negeig)];
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
360 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
361 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
362 end