annotate +scheme/Hypsyst2dCurve.m @ 300:32f3ee81bc37 feature/hypsyst

A commit before i destroy everything
author Ylva Rydin <ylva.rydin@telia.com>
date Tue, 04 Oct 2016 16:25:38 +0200
parents 4d8d6eb0c116
children d9860ebc3148
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
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
8 J, Ji
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
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
21 % Norms in the x and y directions
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
22 Hxii,Hetai % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
23 I_xi,I_eta, I_N, onesN
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
24 e_w, e_e, e_s, e_n
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
25 index_w, index_e,index_s,index_n
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
26 params %parameters for the coeficient matrice
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
27 end
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
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
30 methods
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
31 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
32 default_arg('E', [])
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
33 xilim = {0 1};
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
34 etalim = {0 1};
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
35
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
36 if length(m) == 1
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
37 m = [m m];
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
38 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
39 obj.params = params;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
40 obj.A=A;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
41 obj.B=B;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
42
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
43
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
44 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
45 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
46 obj.E=@(params,x,y,~,~)E(params,x,y);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
47
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
48 m_xi = m(1);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
49 m_eta = m(2);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
50 m_tot=m_xi*m_eta;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
51
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
52 ops_xi = sbp.D2Standard(m_xi,xilim,order);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
53 ops_eta = sbp.D2Standard(m_eta,etalim,order);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
54
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
55 obj.xi = ops_xi.x;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
56 obj.eta = ops_eta.x;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
57
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
58 obj.Xi = kr(obj.xi,ones(m_eta,1));
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
59 obj.Eta = kr(ones(m_xi,1),obj.eta);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
60
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
61 obj.n = length(A(obj.params,0,0));
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
62 obj.onesN=ones(obj.n);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
63
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
64 obj.index_w=1:m_xi;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
65 obj.index_e=(m_tot-m_xi+1):m_tot;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
66 obj.index_s=1:m_xi:(m_tot-m_xi+1);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
67 obj.index_n=(m_xi):m_xi:m_tot;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
68
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
69 I_n = eye(obj.n);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
70 I_xi = speye(m_xi);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
71 obj.I_xi = I_xi;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
72 I_eta = speye(m_eta);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
73 obj.I_eta = I_eta;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
74
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
75 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
76 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
77 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
78 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
79
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
80 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
81 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
82 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
83 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
84
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
85 [X,Y] = ti.map(obj.xi,obj.eta);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
86
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
87 [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
88 [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
89
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
90 obj.X=reshape(X,m_tot,1);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
91 obj.Y=reshape(Y,m_tot,1);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
92 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
93 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
94 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
95 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
96
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
97 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
98 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
99 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
100
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
101 obj.m=m;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
102 obj.h=[ops_xi.h ops_eta.h];
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
103 obj.order=order;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
104 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
105 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
106
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
107 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
108
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
109 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
110
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
111 % 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
112 % 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
113 % 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
114 % 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
115 % 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
116 function [closure, penalty] = boundary_condition(obj,boundary,type,L)
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
117 default_arg('type','char');
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
118 switch type
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
119 case{'c','char'}
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
120 [closure,penalty]=boundary_condition_char(obj,boundary);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
121 case{'general'}
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
122 [closure,penalty]=boundary_condition_general(obj,boundary,L);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
123 otherwise
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
124 error('No such boundary condition')
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
125 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
126 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
127
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
128 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
129 error('An interface function does not exist yet');
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
130 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
131
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
132 function N = size(obj)
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
133 N = obj.m;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
134 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
135
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
136 function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y,x_,y_)
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
137 params=obj.params;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
138
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
139 if isa(mat,'function_handle')
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
140 [rows,cols]=size(mat(params,0,0,0,0));
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
141 x_=kr(x_,obj.onesN);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
142 y_=kr(y_,obj.onesN);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
143 matVec=mat(params,X',Y',x_',y_');
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
144 matVec=sparse(matVec);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
145 side=max(length(X),length(Y));
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
146 else
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
147 matVec=mat;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
148 [rows,cols]=size(matVec);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
149 side=max(length(X),length(Y));
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
150 cols=cols/side;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
151 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
152 ret=kron(ones(rows,cols),speye(side));
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
153
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
154 for ii=1:rows
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
155 for jj=1:cols
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
156 ret((ii-1)*side+1:ii*side,(jj-1)*side+1:jj*side)=diag(matVec(ii,(jj-1)*side+1:jj*side));
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
157 end
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 end
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
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
162 function [closure, penalty]=boundary_condition_char(obj,boundary)
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
163 params=obj.params;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
164 X=obj.X; Y=obj.Y;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
165 xi=obj.xi; eta=obj.eta;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
166 side=max(length(xi),length(eta));
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
167
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
168 switch boundary
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
169 case {'w','W','west'}
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
170 e_=obj.e_w;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
171 mat=obj.Ahat;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
172 boundPos='l';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
173 Hi=obj.Hxii;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
174 [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));
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));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
181 case {'s','S','south'}
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
182 e_=obj.e_s;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
183 mat=obj.Bhat;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
184 boundPos='l';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
185 Hi=obj.Hetai;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
186 [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));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
187 case {'n','N','north'}
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
188 e_=obj.e_n;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
189 mat=obj.Bhat;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
190 boundPos='r';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
191 Hi=obj.Hetai;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
192 [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));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
193 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
194
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
195 pos=signVec(1); zeroval=signVec(2); neg=signVec(3);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
196
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
197 switch boundPos
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
198 case {'l'}
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
199 tau=sparse(obj.n*side,pos);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
200 Vi_plus=Vi(1:pos,:);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
201 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
202 closure=Hi*e_*V*tau*Vi_plus*e_';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
203 penalty=-Hi*e_*V*tau*Vi_plus;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
204 case {'r'}
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
205 tau=sparse(obj.n*side,neg);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
206 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
207 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
208 closure=Hi*e_*V*tau*Vi_minus*e_';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
209 penalty=-Hi*e_*V*tau*Vi_minus;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
210 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
211 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
212
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
213
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
214 function [closure,penalty]=boundary_condition_general(obj,boundary,L)
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
215 params=obj.params;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
216 X=obj.X; Y=obj.Y;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
217 side=max(length(xi),length(eta));
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
218
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
219 switch boundary
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
220 case {'w','W','west'}
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
221 e_=obj.e_w;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
222 mat=obj.Ahat;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
223 boundPos='l';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
224 Hi=obj.Hxii;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
225 [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));
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
226 L=obj.evaluateCoefficientMatrix(L,X(obj.index_w),Y(obj.index_w));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
227 case {'e','E','east'}
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
228 e_=obj.e_e;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
229 mat=obj.Ahat;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
230 boundPos='r';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
231 Hi=obj.Hxii;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
232 [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));
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
233 L=obj.evaluateCoefficientMatrix(L,X(obj.index_e),Y(obj.index_e),[],[]);
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
234 case {'s','S','south'}
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
235 e_=obj.e_s;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
236 mat=obj.Bhat;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
237 boundPos='l';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
238 Hi=obj.Hetai;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
239 [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));
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
240 L=obj.evaluateCoefficientMatrix(L,X(obj.index_s),Y(obj.index_s),[],[]);
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
241 case {'n','N','north'}
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
242 e_=obj.e_n;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
243 mat=obj.Bhat;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
244 boundPos='r';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
245 Hi=obj.Hetai;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
246 [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));
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
247 L=obj.evaluateCoefficientMatrix(L,X(obj.index_n),Y(obj.index_n));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
248 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
249
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
250 pos=signVec(1); zeroval=signVec(2); neg=signVec(3);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
251
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
252 switch boundPos
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
253 case {'l'}
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
254 tau=sparse(obj.n*side,pos);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
255 Vi_plus=Vi(1:pos,:);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
256 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
257 V_plus=V(:,1:pos);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
258 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
259
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
260 tau(1:pos*side,:)=-abs(D(1:pos,1:pos));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
261 R=-inv(L*V_plus)*(L*V_minus);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
262 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
263 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
264 case {'r'}
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
265 tau=sparse(obj.n*side,neg);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
266 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
267 Vi_plus=Vi(1:pos,:);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
268 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
269
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
270 V_plus=V(:,1:pos);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
271 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
272 R=-inv(L*V_minus)*(L*V_plus);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
273 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
274 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
275 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
276 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
277
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
278
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
279 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
280 params=obj.params;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
281 syms xs ys
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
282 if(sum(abs(x_))~=0)
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
283 syms xs_
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
284 else
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
285 xs_=0;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
286 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
287
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
288 if(sum(abs(y_))~=0)
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
289 syms ys_;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
290 else
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
291 ys_=0;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
292 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
293
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
294 [V, D]=eig(mat(params,xs,ys,xs_,ys_));
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
295 syms xs ys xs_ ys_
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
296
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
297 xs=x;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
298 ys=y;
299
4d8d6eb0c116 A commit before I change the boundary treatment
Ylva Rydin <ylva.rydin@telia.com>
parents: 298
diff changeset
299 xs_=x_;
4d8d6eb0c116 A commit before I change the boundary treatment
Ylva Rydin <ylva.rydin@telia.com>
parents: 298
diff changeset
300 ys_=y_;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
301
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
302 side=max(length(x),length(y));
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
303 Dret=zeros(obj.n,side*obj.n);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
304 Vret=zeros(obj.n,side*obj.n);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
305 for ii=1:obj.n
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
306 for jj=1:obj.n
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
307 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
308 Vret(jj,(ii-1)*side+1:side*ii)=eval(V(jj,ii));
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
309 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
310 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
311
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
312 D=sparse(Dret);
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
313 V=sparse(Vret);
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
314 V=obj.evaluateCoefficientMatrix(V,x,y,x_,y_);
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
315 D=obj.evaluateCoefficientMatrix(D,x,y,x_,y_);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
316 DD=diag(D);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
317
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
318 poseig=(DD>0);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
319 zeroeig=(DD==0);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
320 negeig=(DD<0);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
321
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
322 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
323 V=[V(:,poseig) V(:,zeroeig) V(:,negeig)];
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
324 Vi=inv(V);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
325 signVec=[sum(poseig),sum(zeroeig),sum(negeig)];
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
326 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
327
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
328 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
329 end