Mercurial > repos > public > sbplib

annotate +scheme/Hypsyst2dCurve.m @ 774:66eb4a2bbb72 feature/grids

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Remove default scaling of the system. The scaling doens't seem to help actual solutions. One example that fails in the flexural code. With large timesteps the solutions seems to blow up. One particular example is profilePresentation on the tdb_presentation_figures branch with k = 0.0005
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 18 Jul 2018 15:42:52 -0700
parents 9d1fc984f40d
children 459eeb99130f
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
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
4 n % size of system
298
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
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
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
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
13 X_eta, Y_eta
298
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
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
16
298
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
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
19
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
20 H % Discrete norm
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
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
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
25 params % Parameters for the coeficient matrice
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
26 end
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
27
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
28
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
29 methods
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
30 % Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Eu
298
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};
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
35
298
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;
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
42
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
43 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
44 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
45 obj.E=@(params,x,y,~,~)E(params,x,y);
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
46
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
47 m_xi = m(1);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
48 m_eta = m(2);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
49 m_tot=m_xi*m_eta;
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
50
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
51 ops_xi = sbp.D2Standard(m_xi,xilim,order);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
52 ops_eta = sbp.D2Standard(m_eta,etalim,order);
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
53
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
54 obj.xi = ops_xi.x;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
55 obj.eta = ops_eta.x;
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
56
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
57 obj.Xi = kr(obj.xi,ones(m_eta,1));
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
58 obj.Eta = kr(ones(m_xi,1),obj.eta);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
59
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
60 obj.n = length(A(obj.params,0,0));
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
61 obj.onesN=ones(obj.n);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
62
301
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
63 obj.index_w=1:m_eta;
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
64 obj.index_e=(m_tot-m_e
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
65
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
66 metric_termsta+1):m_tot;
301
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 300
diff changeset
67 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
68 obj.index_n=(m_eta):m_eta:m_tot;
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
69
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
70 I_n = eye(obj.n);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
71 I_xi = speye(m_xi);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
72 obj.I_xi = I_xi;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
73 I_eta = speye(m_eta);
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
74 obj.I_eta = I_eta;
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
75
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
76 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
77 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
78 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
79 obj.Hetai = kr(I_n, I_xi, ops_eta.HI);
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
80
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
81 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
82 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
83 obj.e_s = kr(I_n, I_xi, ops_eta.e_l);
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
84 obj.e_n = kr(I_n, I_xi,
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
85
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
86 metric_termsops_eta.e_r);
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
87
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
88 [X,Y] = ti.map(obj.xi,obj.eta);
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
89
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
90 [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
91 [y_xi,y_eta] = gridDerivatives(Y,ops_xi.D1, ops_eta.D1);
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
92
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
93 obj.X = reshape(X,m_tot,1);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
94 obj.Y = reshape(Y,m_tot,1);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
95 obj.X_xi = reshape(x_xi,m_tot,1);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
96 obj.Y_xi = reshape(y_xi,m_tot,1);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
97 obj.X_eta = reshape(x_eta,m_tot,1);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
98 obj.Y_eta = reshape(y_eta,m_tot,1);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
99
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
100 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
101 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
102 E_evaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,[],[]);
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
103
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
104 obj.m = m;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
105 obj.h = [ops_xi.h ops_eta.h];
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
106 obj.order = order;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
107 obj.J = obj.X_xi.*obj.Y_eta - obj.X_eta.*obj.Y_xi;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
108 obj.Ji = kr(I_n,spdiags(1./obj.J,0,m_tot,m_tot));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
109
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
110 obj.D = obj.Ji*(-Ahat_evaluated*D1_xi-Bhat_evaluated*D1_eta)-E_evaluated;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
111
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
112 end
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
113
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
114 % 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
115 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
116 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w',General boundary conditions'n','s'.
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
117 % 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
118 % 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
119 function [closure, penalty] = boundary_condition(obj,boundary,type,L)
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
120 default_arg('type','char');
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
121 switch type
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
122 case{'c','char'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
123 [closure,penalty] = boundary_condition_char(obj,boundary);
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
124 case{'general'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
125 [closure,penalty] = boundary_condition_general(obj,boundary,L);
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
126 otherwise
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
127 error('No such boundary condition')
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 end
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
130
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
131 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundaryGeneral boundary conditions)
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
132 error('An interface function does not exist yet');
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
133 end
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
134
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
135 function N = size(obj)
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
136 N = obj.m;
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
137 end
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
138
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
139 function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y,x_,y_)
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
140 params = obj.params;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
141
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
142 if isa(mat,'function_handle')
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
143 [rows,cols] = size(mat(params,0,0,0,0));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
144 x_ = kr(obj.onesN,x_);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
145 y_ = kr(obj.onesN,y_);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
146 matVec = mat(params,X',Y',x_',y_');
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
147 matVec = sparse(matVec);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
148 side = max(length(X),length(Y));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
149 else
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
150 matVec = mat;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
151 [rows,cols] = size(matVec);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
152 side = max(length(X),length(Y));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
153 cols = cols/side;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
154 end
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
155
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
156 ret = cell(rows,cols);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
157 for ii = 1:rows
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
158 for jj = 1:cols
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
159 ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
160 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
161 end
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
162 ret = cell2mat(ret);
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
163 end
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
164
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
165 %Characteristic boundary conditions
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
166 function [closure, penalty] = boundary_condition_char(obj,boundary)
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
167 params = obj.params;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
168 X = obj.X;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
169 Y = obj.Y;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
170 xi = obj.xi;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
171 eta = obj.eta;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
172
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
173 switch boundary
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
174 case {'w','W','west'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
175 e_ = obj.e_w;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
176 mat = obj.Ahat;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
177 boundPos = 'l';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
178 Hi = obj.Hxii;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
179 [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));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
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 {'e','E','east'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
182 e_ = obj.e_e;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
183 mat = obj.Ahat;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
184 boundPos = 'r';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
185 Hi = obj.Hxii;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
186 [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));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
187 side = max(length(eta));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
188 case {'s','S','south'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
189 e_ = obj.e_s;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
190 mat = obj.Bhat;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
191 boundPos = 'l';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
192 Hi = obj.Hetai;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
193 [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));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
194 side = max(length(xi));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
195 case {'n','N','north'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
196 e_ = obj.e_n;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
197 mat = obj.Bhat;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
198 boundPos = 'r';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
199 Hi = obj.Hetai;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
200 [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));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
201 side = max(length(xi));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
202 end
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
203
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
204 pos = signVec(1);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
205 zeroval = signVec(2);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
206 neg = signVec(3);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
207
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
208 switch boundPos
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
209 case {'l'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
210 tau = sparse(obj.n*side,pos);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
211 Vi_plus = Vi(1:pos,:);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
212 tau(1:pos,:) = -abs(D(1:pos,1:pos));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
213 closure = Hi*e_*V*tau*Vi_plus*e_';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
214 penalty = -Hi*e_*V*tau*Vi_plus;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
215 case {'r'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
216 tau = sparse(obj.n*side,neg);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
217 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
218 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
219 closure = Hi*e_*V*tau*Vi_minus*e_';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
220 penalty = -Hi*e_*V*tau*Vi_minus;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
221 end
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
222 end
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
223
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
224
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
225 % General boundary condition in the form Lu=g(x)
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
226 function [closure,penalty] = boundary_condition_general(obj,boundary,L)
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
227 params = obj.params;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
228 X = obj.X;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
229 Y = obj.Y;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
230 xi = obj.xi;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
231 eta = obj.eta;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
232
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
233 switch boundary
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
234 case {'w','W','west'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
235 e_ = obj.e_w;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
236 mat = obj.Ahat;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
237 boundPos = 'l';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
238 Hi = obj.Hxii;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
239 [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));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
240
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
241 Ji_vec = diag(obj.Ji);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
242 Ji = diag(Ji_vec(obj.index_w));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
243 xi_x = Ji*obj.Y_eta(obj.index_w);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
244 xi_y = -Ji*obj.X_eta(obj.index_w);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
245 L = obj.evaluateCoefficientMatrix(L,xi_x,xi_y,[],[]);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
246 side = max(length(eta));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
247 case {'e','E','east'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
248 e_ = obj.e_e;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
249 mat = obj.Ahat;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
250 boundPos = 'r';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
251 Hi = obj.Hxii;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
252 [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
253
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
254 Ji_vec = diag(obj.Ji);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
255 Ji = diag(Ji_vec(obj.index_e));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
256 xi_x = Ji*obj.Y_eta(obj.index_e);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
257 xi_y = -Ji*obj.X_eta(obj.index_e);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
258 L = obj.evaluateCoefficientMatrix(L,-xi_x,-xi_y,[],[]);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
259 side = max(length(eta));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
260 case {'s','S','south'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
261 e_ = obj.e_s;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
262 mat = obj.Bhat;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
263 boundPos = 'l';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
264 Hi = obj.Hetai;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
265 [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
266
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
267 Ji_vec = diag(obj.Ji);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
268 Ji = diag(Ji_vec(obj.index_s));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
269 eta_x = Ji*obj.Y_xi(obj.index_s);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
270 eta_y = -Ji*obj.X_xi(obj.index_s);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
271 L = obj.evaluateCoefficientMatrix(L,eta_x,eta_y,[],[]);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
272 side = max(length(xi));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
273 case {'n','N','north'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
274 e_ = obj.e_n;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
275 mat = obj.Bhat;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
276 boundPos = 'r';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
277 Hi = obj.Hetai;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
278 [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
279
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
280 Ji_vec = diag(obj.Ji);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
281 Ji = diag(Ji_vec(obj.index_n));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
282 eta_x = Ji*obj.Y_xi(obj.index_n);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
283 eta_y = -Ji*obj.X_xi(obj.index_n);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
284 L = obj.evaluateCoefficientMatrix(L,-eta_x,-eta_y,[],[]);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
285 side = max(length(xi));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
286 end
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
287
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
288 pos = signVec(1);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
289 zeroval = signVec(2);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
290 neg = signVec(3);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
291
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
292 switch boundPos
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
293 case {'l'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
294 tau = sparse(obj.n*side,pos);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
295 Vi_plus = Vi(1:pos,:);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
296 Vi_minus = Vi(pos+1:obj.n*side,:);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
297 V_plus = V(:,1:pos);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
298 V_minus = V(:,(pos)+1:obj.n*side);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
299
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
300 tau(1:pos,:) = -abs(D(1:pos,1:pos));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
301 R = -inv(L*V_plus)*(L*V_minus);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
302 closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
303 penalty = -Hi*e_*V*tau*inv(L*V_plus)*L;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
304 case {'r'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
305 tau = sparse(obj.n*side,neg);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
306 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
307 Vi_plus = Vi(1:pos,:);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
308 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
309
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
310 V_plus = V(:,1:pos);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
311 V_minus = V(:,(pos+zeroval)+1:obj.n*side);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
312 R = -inv(L*V_minus)*(L*V_plus);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
313 closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
314 penalty = -Hi*e_*V*tau*inv(L*V_minus)*L;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
315 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
316 end
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
317
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
318 % Function that diagonalizes a symbolic matrix A as A=V*D*Vi
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
319 % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
320 % [d+ ]
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
321 % D = [ d0 ]
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
322 % [ d-]
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
323 % signVec is a vector specifying the number of possitive, zero and negative eigenvalues of D
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
324 function [V,Vi, D,signVec] = matrixDiag(obj,mat,x,y,x_,y_)
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
325 params = obj.params;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
326 syms xs ys
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
327 if(sum(abs(x_)) ~= 0)
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
328 syms xs_
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
329 else
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
330 xs_ = 0;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
331 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
332
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
333 if(sum(abs(y_))~= 0)
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
334 syms ys_;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
335 else
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
336 ys_ = 0;
298
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
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
339 [V, D] = eig(mat(params,xs,ys,xs_,ys_));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
340 Vi = inv(V);
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
341 syms xs ys xs_ ys_
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
342
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
343 xs = x;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
344 ys = y;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
345 xs_ = x_;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
346 ys_ = y_;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
347
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
348 side = max(length(x),length(y));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
349 Dret = zeros(obj.n,side*obj.n);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
350 Vret = zeros(obj.n,side*obj.n);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
351 Viret = zeros(obj.n,side*obj.n);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
352 for ii = 1:obj.n
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
353 for jj = 1:obj.n
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
354 Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
355 Vret(jj,(ii-1)*side+1:side*ii) = eval(V(jj,ii));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
356 Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
357 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
358 end
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
359
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
360 D = sparse(Dret);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
361 V = sparse(Vret);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
362 Vi = sparse(Viret);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
363 V = obj.evaluateCoefficientMatrix(V,x,y,x_,y_);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
364 D = obj.evaluateCoefficientMatrix(D,x,y,x_,y_);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
365 Vi = obj.evaluateCoefficientMatrix(Vi,x,y,x_,y_);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
366 DD = diag(D);
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
367
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
368 poseig = (DD>0);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
369 zeroeig = (DD==0);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
370 negeig = (DD<0);
300
32f3ee81bc37 A commit before i destroy everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 299
diff changeset
371
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
372 D = diag([DD(poseig); DD(zeroeig); DD(negeig)]);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
373 V = [V(:,poseig) V(:,zeroeig) V(:,negeig)];
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
374 Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
375 signVec = [sum(poseig),sum(zeroeig),sum(negeig)];
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
376 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
377 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
378 end