Mercurial > repos > public > sbplib

annotate +scheme/Hypsyst2dCurve.m @ 301:d9860ebc3148 feature/hypsyst

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
HypsystCurve2D Seems to work (Converges with MMS)
author Ylva Rydin <ylva.rydin@telia.com>
date Wed, 05 Oct 2016 17:36:34 +0200
parents 32f3ee81bc37
children 5d5652fe826a
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
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
150 ret=kron(ones(rows,cols),speye(side));
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
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
154 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
155 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
156 end
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
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 function [closure, penalty]=boundary_condition_char(obj,boundary)
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
161 params=obj.params;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
162 X=obj.X; Y=obj.Y;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
163 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
164
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
165
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
166 switch boundary
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
167 case {'w','W','west'}
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
168 e_=obj.e_w;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
169 mat=obj.Ahat;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
170 boundPos='l';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
171 Hi=obj.Hxii;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
172 [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
173 side=max(length(eta));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
174 case {'e','E','east'}
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
175 e_=obj.e_e;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
176 mat=obj.Ahat;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
177 boundPos='r';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
178 Hi=obj.Hxii;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
179 [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
180 side=max(length(eta));
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));
301
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
187 side=max(length(xi));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
188 case {'n','N','north'}
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
189 e_=obj.e_n;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
190 mat=obj.Bhat;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
191 boundPos='r';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
192 Hi=obj.Hetai;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
193 [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
194 side=max(length(xi));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
195 end
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 pos=signVec(1); zeroval=signVec(2); neg=signVec(3);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
198
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
199 switch boundPos
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
200 case {'l'}
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
201 tau=sparse(obj.n*side,pos);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
202 Vi_plus=Vi(1:pos,:);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
203 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
204 closure=Hi*e_*V*tau*Vi_plus*e_';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
205 penalty=-Hi*e_*V*tau*Vi_plus;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
206 case {'r'}
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
207 tau=sparse(obj.n*side,neg);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
208 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
209 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
210 closure=Hi*e_*V*tau*Vi_minus*e_';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
211 penalty=-Hi*e_*V*tau*Vi_minus;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
212 end
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
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 function [closure,penalty]=boundary_condition_general(obj,boundary,L)
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
217 params=obj.params;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
218 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
219 xi=obj.xi; eta=obj.eta;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
220
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
221 switch boundary
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
222 case {'w','W','west'}
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
223 e_=obj.e_w;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
224 mat=obj.Ahat;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
225 boundPos='l';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
226 Hi=obj.Hxii;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
227 [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
228
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
229 Ji_vec=diag(obj.Ji);
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
230 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
231 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
232 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
233 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
234 side=max(length(eta));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
235 case {'e','E','east'}
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
236 e_=obj.e_e;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
237 mat=obj.Ahat;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
238 boundPos='r';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
239 Hi=obj.Hxii;
301
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
240 [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
241
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
242 Ji_vec=diag(obj.Ji);
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
243 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
244 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
245 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
246 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
247 side=max(length(eta));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
248 case {'s','S','south'}
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
249 e_=obj.e_s;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
250 mat=obj.Bhat;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
251 boundPos='l';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
252 Hi=obj.Hetai;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
253 [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
254
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
255 Ji_vec=diag(obj.Ji);
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
256 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
257 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
258 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
259 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
260 side=max(length(xi));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
261 case {'n','N','north'}
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
262 e_=obj.e_n;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
263 mat=obj.Bhat;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
264 boundPos='r';
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
265 Hi=obj.Hetai;
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
266 [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
267
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
268 Ji_vec=diag(obj.Ji);
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
269 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
270 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
271 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
272 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
273
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
274 side=max(length(xi));
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
275
298
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 pos=signVec(1); zeroval=signVec(2); neg=signVec(3);
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 switch boundPos
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
281 case {'l'}
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
282 tau=sparse(obj.n*side,pos);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
283 Vi_plus=Vi(1:pos,:);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
284 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
285 V_plus=V(:,1:pos);
301
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
286 V_minus=V(:,(pos)+1:obj.n*side);
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
287
301
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
288 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
289 R=-inv(L*V_plus)*(L*V_minus);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
290 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
291 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
292 case {'r'}
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
293 tau=sparse(obj.n*side,neg);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
294 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
295 Vi_plus=Vi(1:pos,:);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
296 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
297
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
298 V_plus=V(:,1:pos);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
299 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
300 R=-inv(L*V_minus)*(L*V_plus);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
301 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
302 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
303 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
304 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
305
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
306
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
307 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
308 params=obj.params;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
309 syms xs ys
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
310 if(sum(abs(x_))~=0)
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
311 syms xs_
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
312 else
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
313 xs_=0;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
314 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
315
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
316 if(sum(abs(y_))~=0)
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
317 syms ys_;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
318 else
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
319 ys_=0;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
320 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
321
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
322 [V, D]=eig(mat(params,xs,ys,xs_,ys_));
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
323 syms xs ys xs_ ys_
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
324
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
325 xs=x;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
326 ys=y;
299
4d8d6eb0c116 A commit before I change the boundary treatment
Ylva Rydin <ylva.rydin@telia.com>
parents: 298
diff changeset
327 xs_=x_;
4d8d6eb0c116 A commit before I change the boundary treatment
Ylva Rydin <ylva.rydin@telia.com>
parents: 298
diff changeset
328 ys_=y_;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
329
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
330 side=max(length(x),length(y));
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
331 Dret=zeros(obj.n,side*obj.n);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
332 Vret=zeros(obj.n,side*obj.n);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
333 for ii=1:obj.n
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
334 for jj=1:obj.n
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
335 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
336 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
337 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
338 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
339
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
340 D=sparse(Dret);
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
341 V=sparse(Vret);
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
342 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
343 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
344 DD=diag(D);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
345
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
346 poseig=(DD>0);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
347 zeroeig=(DD==0);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
348 negeig=(DD<0);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
349
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
350 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
351 V=[V(:,poseig) V(:,zeroeig) V(:,negeig)];
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
352 Vi=inv(V);
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
353 signVec=[sum(poseig),sum(zeroeig),sum(negeig)];
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
354 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
355
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
356 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
357 end