Mercurial > repos > public > sbplib

annotate +scheme/Hypsyst2dCurve.m @ 1198:2924b3a9b921 feature/d2_compatible

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add OpSet for fully compatible D2Variable, created from regular D2Variable by replacing d1 by first row of D1. Formal reduction by one order of accuracy at the boundary point.
author Martin Almquist <malmquist@stanford.edu>
date Fri, 16 Aug 2019 14:30:28 -0700
parents 8d73fcdb07a5
children
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
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
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
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
11
298
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
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
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
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
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
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
27
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
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};
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
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;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
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);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
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;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
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);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
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;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
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);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
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);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
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;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
64 obj.index_e=(m_tot-m_e
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
65
369
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;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
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;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
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);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
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);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
84 obj.e_n = kr(I_n, I_xi,
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
85
369
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);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
87
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
88 [X,Y] = ti.map(obj.xi,obj.eta);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
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);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
92
369
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);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
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,[],[]);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
103
369
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));
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
109
369
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;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
111
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
112 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
113
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
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
130
946
706d1c2b4199 Raname opts to type in a bunch of interface methods
Jonatan Werpers <jonatan@werpers.com>
parents: 910
diff changeset
131 function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type)
706d1c2b4199 Raname opts to type in a bunch of interface methods
Jonatan Werpers <jonatan@werpers.com>
parents: 910
diff changeset
132 error('Not implemented');
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
133 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
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
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
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;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
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
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
155
369
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
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
164
369
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;
997
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
172 e_ = obj.getBoundaryOperator('e', boundary);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
173
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
174 switch boundary
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
175 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
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 mat = obj.Ahat;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
183 boundPos = 'r';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
184 Hi = obj.Hxii;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
185 [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
186 side = max(length(eta));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
187 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
188 mat = obj.Bhat;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
189 boundPos = 'l';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
190 Hi = obj.Hetai;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
191 [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
192 side = max(length(xi));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
193 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
194 mat = obj.Bhat;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
195 boundPos = 'r';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
196 Hi = obj.Hetai;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
197 [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
198 side = max(length(xi));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
199 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
200
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
201 pos = signVec(1);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
202 zeroval = signVec(2);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
203 neg = signVec(3);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
204
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
205 switch boundPos
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
206 case {'l'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
207 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
208 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
209 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
210 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
211 penalty = -Hi*e_*V*tau*Vi_plus;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
212 case {'r'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
213 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
214 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
215 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
216 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
217 penalty = -Hi*e_*V*tau*Vi_minus;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
218 end
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
219 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
220
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
221
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
222 % 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
223 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
224 params = obj.params;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
225 X = obj.X;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
226 Y = obj.Y;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
227 xi = obj.xi;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
228 eta = obj.eta;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
229
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
230 switch boundary
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
231 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
232 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
233 mat = obj.Ahat;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
234 boundPos = 'l';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
235 Hi = obj.Hxii;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
236 [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));
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
237
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
238 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
239 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
240 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
241 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
242 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
243 side = max(length(eta));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
244 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
245 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
246 mat = obj.Ahat;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
247 boundPos = 'r';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
248 Hi = obj.Hxii;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
249 [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));
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
250
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
251 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
252 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
253 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
254 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
255 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
256 side = max(length(eta));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
257 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
258 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
259 mat = obj.Bhat;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
260 boundPos = 'l';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
261 Hi = obj.Hetai;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
262 [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));
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
263
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
264 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
265 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
266 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
267 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
268 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
269 side = max(length(xi));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
270 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
271 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
272 mat = obj.Bhat;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
273 boundPos = 'r';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
274 Hi = obj.Hetai;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
275 [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));
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
276
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
277 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
278 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
279 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
280 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
281 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
282 side = max(length(xi));
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
283 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
284
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
285 pos = signVec(1);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
286 zeroval = signVec(2);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
287 neg = signVec(3);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
288
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
289 switch boundPos
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
290 case {'l'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
291 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
292 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
293 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
294 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
295 V_minus = V(:,(pos)+1:obj.n*side);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
296
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
297 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
298 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
299 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
300 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
301 case {'r'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
302 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
303 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
304 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
305 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
306
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
307 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
308 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
309 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
310 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
311 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
312 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
313 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
314
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
315 % 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
316 % 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
317 % [d+ ]
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
318 % D = [ d0 ]
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
319 % [ d-]
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
320 % 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
321 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
322 params = obj.params;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
323 syms xs ys
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
324 if(sum(abs(x_)) ~= 0)
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
325 syms xs_
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
326 else
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
327 xs_ = 0;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
328 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
329
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
330 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
331 syms ys_;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
332 else
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
333 ys_ = 0;
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
334 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
335
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
336 [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
337 Vi = inv(V);
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
338 syms xs ys xs_ ys_
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
339
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
340 xs = x;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
341 ys = y;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
342 xs_ = x_;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
343 ys_ = y_;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
344
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
345 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
346 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
347 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
348 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
349 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
350 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
351 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
352 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
353 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
354 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
355 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
356
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
357 D = sparse(Dret);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
358 V = sparse(Vret);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
359 Vi = sparse(Viret);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
360 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
361 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
362 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
363 DD = diag(D);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
364
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
365 poseig = (DD>0);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
366 zeroeig = (DD==0);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
367 negeig = (DD<0);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
368
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
369 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
370 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
371 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
372 signVec = [sum(poseig),sum(zeroeig),sum(negeig)];
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
373 end
997
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
374
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
375 % Returns the boundary operator op for the boundary specified by the string boundary.
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
376 % op -- string or a cell array of strings
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
377 % boundary -- string
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
378 function varargout = getBoundaryOperator(obj, op, boundary)
1042
8d73fcdb07a5 Add asserts to boundary identifier inputs
Jonatan Werpers <jonatan@werpers.com>
parents: 997
diff changeset
379 assertIsMember(boundary, {'w', 'e', 's', 'n'})
997
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
380
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
381 if ~iscell(op)
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
382 op = {op};
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
383 end
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
384
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
385 for i = 1:numel(op)
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
386 switch op{i}
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
387 case 'e'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
388 switch boundary
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
389 case 'w'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
390 e = obj.e_w;
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
391 case 'e'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
392 e = obj.e_e;
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
393 case 's'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
394 e = obj.e_s;
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
395 case 'n'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
396 e = obj.e_n;
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
397 end
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
398 varargout{i} = e;
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
399 end
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
400 end
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
401 end
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
402
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
403 % Returns square boundary quadrature matrix, of dimension
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
404 % corresponding to the number of boundary points
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
405 %
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
406 % boundary -- string
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
407 function H_b = getBoundaryQuadrature(obj, boundary)
1042
8d73fcdb07a5 Add asserts to boundary identifier inputs
Jonatan Werpers <jonatan@werpers.com>
parents: 997
diff changeset
408 assertIsMember(boundary, {'w', 'e', 's', 'n'})
997
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
409
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
410 e = obj.getBoundaryOperator('e', boundary);
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
411
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
412 switch boundary
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
413 case 'w'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
414 H_b = inv(e'*obj.Hetai*e);
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
415 case 'e'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
416 H_b = inv(e'*obj.Hetai*e);
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
417 case 's'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
418 H_b = inv(e'*obj.Hxii*e);
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
419 case 'n'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
420 H_b = inv(e'*obj.Hxii*e);
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
421 end
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
422 end
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
423
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
424
298
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
425 end
861972361f75 The curvelinear works for quads
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
426 end