Mercurial > repos > public > sbplib

annotate +scheme/Heat2dVariable.m @ 1303:49e3870335ef feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make the hollow scheme generation more efficient by introducing the D2VariableHollow opSet
author Martin Almquist <malmquist@stanford.edu>
date Sat, 11 Jul 2020 06:54:15 -0700
parents 21394c78c72e
children 78db023a7fe3
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
688
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
1 classdef Heat2dVariable < scheme.Scheme
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
2
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
3 % Discretizes the Laplacian with variable coefficent,
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
4 % In the Heat equation way (i.e., the discretization matrix is not necessarily
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
5 % symmetric)
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
6 % u_t = div * (kappa * grad u )
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
7 % opSet should be cell array of opSets, one per dimension. This
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
8 % is useful if we have periodic BC in one direction.
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
9
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
10 properties
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
11 m % Number of points in each direction, possibly a vector
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
12 h % Grid spacing
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
13
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
14 grid
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
15 dim
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
16
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
17 order % Order of accuracy for the approximation
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
18
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
19 % Diagonal matrix for variable coefficients
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
20 KAPPA % Variable coefficient
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
21
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
22 D % Total operator
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
23 D1 % First derivatives
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
24
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
25 % Second derivatives
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
26 D2_kappa
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
27
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
28 H, Hi % Inner products
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
29 e_l, e_r
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
30 d1_l, d1_r % Normal derivatives at the boundary
725
e9e15d64f9f9 Implement symmetric Dirichlet BC with variable tuning parameter in Heat2D scheme.
Martin Almquist <malmquist@stanford.edu>
parents: 724
diff changeset
31 alpha % Vector of borrowing constants
688
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
32
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
33 H_boundary % Boundary inner products
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
34
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
35 end
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
36
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
37 methods
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
38
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
39 function obj = Heat2dVariable(g ,order, kappa_fun, opSet)
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
40 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
41 default_arg('kappa_fun', @(x,y) 0*x+1);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
42 dim = 2;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
43
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
44 assert(isa(g, 'grid.Cartesian'))
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
45
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
46 kappa = grid.evalOn(g, kappa_fun);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
47 m = g.size();
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
48 m_tot = g.N();
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
49
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
50 h = g.scaling();
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
51 lim = g.lim;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
52
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
53 % 1D operators
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
54 ops = cell(dim,1);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
55 for i = 1:dim
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
56 ops{i} = opSet{i}(m(i), lim{i}, order);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
57 end
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
58
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
59 I = cell(dim,1);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
60 D1 = cell(dim,1);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
61 D2 = cell(dim,1);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
62 H = cell(dim,1);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
63 Hi = cell(dim,1);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
64 e_l = cell(dim,1);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
65 e_r = cell(dim,1);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
66 d1_l = cell(dim,1);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
67 d1_r = cell(dim,1);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
68
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
69 for i = 1:dim
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
70 I{i} = speye(m(i));
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
71 D1{i} = ops{i}.D1;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
72 D2{i} = ops{i}.D2;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
73 H{i} = ops{i}.H;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
74 Hi{i} = ops{i}.HI;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
75 e_l{i} = ops{i}.e_l;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
76 e_r{i} = ops{i}.e_r;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
77 d1_l{i} = ops{i}.d1_l;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
78 d1_r{i} = ops{i}.d1_r;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
79 end
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
80
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
81 %====== Assemble full operators ========
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
82 KAPPA = spdiag(kappa);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
83 obj.KAPPA = KAPPA;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
84
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
85 obj.D1 = cell(dim,1);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
86 obj.D2_kappa = cell(dim,1);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
87 obj.e_l = cell(dim,1);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
88 obj.e_r = cell(dim,1);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
89 obj.d1_l = cell(dim,1);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
90 obj.d1_r = cell(dim,1);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
91
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
92 % D1
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
93 obj.D1{1} = kron(D1{1},I{2});
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
94 obj.D1{2} = kron(I{1},D1{2});
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
95
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
96 % Boundary operators
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
97 obj.e_l{1} = kron(e_l{1},I{2});
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
98 obj.e_l{2} = kron(I{1},e_l{2});
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
99 obj.e_r{1} = kron(e_r{1},I{2});
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
100 obj.e_r{2} = kron(I{1},e_r{2});
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
101
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
102 obj.d1_l{1} = kron(d1_l{1},I{2});
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
103 obj.d1_l{2} = kron(I{1},d1_l{2});
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
104 obj.d1_r{1} = kron(d1_r{1},I{2});
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
105 obj.d1_r{2} = kron(I{1},d1_r{2});
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
106
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
107 % D2
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
108 for i = 1:dim
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
109 obj.D2_kappa{i} = sparse(m_tot);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
110 end
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
111 ind = grid.funcToMatrix(g, 1:m_tot);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
112
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
113 for i = 1:m(2)
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
114 D_kappa = D2{1}(kappa(ind(:,i)));
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
115 p = ind(:,i);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
116 obj.D2_kappa{1}(p,p) = D_kappa;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
117 end
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
118
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
119 for i = 1:m(1)
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
120 D_kappa = D2{2}(kappa(ind(i,:)));
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
121 p = ind(i,:);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
122 obj.D2_kappa{2}(p,p) = D_kappa;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
123 end
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
124
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
125 % Quadratures
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
126 obj.H = kron(H{1},H{2});
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
127 obj.Hi = inv(obj.H);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
128 obj.H_boundary = cell(dim,1);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
129 obj.H_boundary{1} = H{2};
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
130 obj.H_boundary{2} = H{1};
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
131
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
132 % Differentiation matrix D (without SAT)
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
133 D2_kappa = obj.D2_kappa;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
134 D1 = obj.D1;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
135 D = sparse(m_tot,m_tot);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
136 for i = 1:dim
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
137 D = D + D2_kappa{i};
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
138 end
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
139 obj.D = D;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
140 %=========================================%
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
141
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
142 % Misc.
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
143 obj.m = m;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
144 obj.h = h;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
145 obj.order = order;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
146 obj.grid = g;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
147 obj.dim = dim;
725
e9e15d64f9f9 Implement symmetric Dirichlet BC with variable tuning parameter in Heat2D scheme.
Martin Almquist <malmquist@stanford.edu>
parents: 724
diff changeset
148 obj.alpha = [ops{1}.borrowing.M.d1, ops{2}.borrowing.M.d1];
688
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
149
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
150 end
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
151
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
152
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
153 % Closure functions return the operators applied to the own domain to close the boundary
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
154 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
155 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
156 % type is a string specifying the type of boundary condition.
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
157 % data is a function returning the data that should be applied at the boundary.
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
158 % neighbour_scheme is an instance of Scheme that should be interfaced to.
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
159 % neighbour_boundary is a string specifying which boundary to interface to.
725
e9e15d64f9f9 Implement symmetric Dirichlet BC with variable tuning parameter in Heat2D scheme.
Martin Almquist <malmquist@stanford.edu>
parents: 724
diff changeset
160 function [closure, penalty] = boundary_condition(obj, boundary, type, symmetric, tuning)
688
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
161 default_arg('type','Neumann');
725
e9e15d64f9f9 Implement symmetric Dirichlet BC with variable tuning parameter in Heat2D scheme.
Martin Almquist <malmquist@stanford.edu>
parents: 724
diff changeset
162 default_arg('symmetric', false);
e9e15d64f9f9 Implement symmetric Dirichlet BC with variable tuning parameter in Heat2D scheme.
Martin Almquist <malmquist@stanford.edu>
parents: 724
diff changeset
163 default_arg('tuning',1.2);
688
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
164
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
165 % j is the coordinate direction of the boundary
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
166 % nj: outward unit normal component.
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
167 % nj = -1 for west, south, bottom boundaries
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
168 % nj = 1 for east, north, top boundaries
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
169 [j, nj] = obj.get_boundary_number(boundary);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
170 switch nj
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
171 case 1
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
172 e = obj.e_r;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
173 d = obj.d1_r;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
174 case -1
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
175 e = obj.e_l;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
176 d = obj.d1_l;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
177 end
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
178
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
179 Hi = obj.Hi;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
180 H_gamma = obj.H_boundary{j};
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
181 KAPPA = obj.KAPPA;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
182 kappa_gamma = e{j}'*KAPPA*e{j};
725
e9e15d64f9f9 Implement symmetric Dirichlet BC with variable tuning parameter in Heat2D scheme.
Martin Almquist <malmquist@stanford.edu>
parents: 724
diff changeset
183 h = obj.h(j);
e9e15d64f9f9 Implement symmetric Dirichlet BC with variable tuning parameter in Heat2D scheme.
Martin Almquist <malmquist@stanford.edu>
parents: 724
diff changeset
184 alpha = h*obj.alpha(j);
688
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
185
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
186 switch type
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
187
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
188 % Dirichlet boundary condition
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
189 case {'D','d','dirichlet','Dirichlet'}
725
e9e15d64f9f9 Implement symmetric Dirichlet BC with variable tuning parameter in Heat2D scheme.
Martin Almquist <malmquist@stanford.edu>
parents: 724
diff changeset
190
e9e15d64f9f9 Implement symmetric Dirichlet BC with variable tuning parameter in Heat2D scheme.
Martin Almquist <malmquist@stanford.edu>
parents: 724
diff changeset
191 if ~symmetric
688
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
192 closure = -nj*Hi*d{j}*kappa_gamma*H_gamma*(e{j}' );
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
193 penalty = nj*Hi*d{j}*kappa_gamma*H_gamma;
725
e9e15d64f9f9 Implement symmetric Dirichlet BC with variable tuning parameter in Heat2D scheme.
Martin Almquist <malmquist@stanford.edu>
parents: 724
diff changeset
194 else
e9e15d64f9f9 Implement symmetric Dirichlet BC with variable tuning parameter in Heat2D scheme.
Martin Almquist <malmquist@stanford.edu>
parents: 724
diff changeset
195 closure = nj*Hi*d{j}*kappa_gamma*H_gamma*(e{j}' )...
e9e15d64f9f9 Implement symmetric Dirichlet BC with variable tuning parameter in Heat2D scheme.
Martin Almquist <malmquist@stanford.edu>
parents: 724
diff changeset
196 -tuning*2/alpha*Hi*e{j}*kappa_gamma*H_gamma*(e{j}' ) ;
e9e15d64f9f9 Implement symmetric Dirichlet BC with variable tuning parameter in Heat2D scheme.
Martin Almquist <malmquist@stanford.edu>
parents: 724
diff changeset
197 penalty = -nj*Hi*d{j}*kappa_gamma*H_gamma ...
e9e15d64f9f9 Implement symmetric Dirichlet BC with variable tuning parameter in Heat2D scheme.
Martin Almquist <malmquist@stanford.edu>
parents: 724
diff changeset
198 +tuning*2/alpha*Hi*e{j}*kappa_gamma*H_gamma;
e9e15d64f9f9 Implement symmetric Dirichlet BC with variable tuning parameter in Heat2D scheme.
Martin Almquist <malmquist@stanford.edu>
parents: 724
diff changeset
199 end
688
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
200
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
201 % Free boundary condition
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
202 case {'N','n','neumann','Neumann'}
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
203 closure = -nj*Hi*e{j}*kappa_gamma*H_gamma*(d{j}' );
724
a9e8c9d71307 Modify penalty for Neumann in Head2d so that data is for normal derivative and not u_x or u_y.
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
204 penalty = Hi*e{j}*kappa_gamma*H_gamma;
a9e8c9d71307 Modify penalty for Neumann in Head2d so that data is for normal derivative and not u_x or u_y.
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
205 % penalty is for normal derivative and not for derivative, hence the sign.
688
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
206
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
207 % Unknown boundary condition
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
208 otherwise
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
209 error('No such boundary condition: type = %s',type);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
210 end
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
211 end
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
212
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
213 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
214 % u denotes the solution in the own domain
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
215 % v denotes the solution in the neighbour domain
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
216 error('Interface not implemented');
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
217 end
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
218
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
219 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
220 function [j, nj] = get_boundary_number(obj, boundary)
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
221
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
222 switch boundary
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
223 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
224 j = 1;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
225 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
226 j = 2;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
227 otherwise
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
228 error('No such boundary: boundary = %s',boundary);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
229 end
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
230
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
231 switch boundary
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
232 case {'w','W','west','West','s','S','south','South'}
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
233 nj = -1;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
234 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
235 nj = 1;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
236 end
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
237 end
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
238
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
239 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
240 function [return_op] = get_boundary_operator(obj, op, boundary)
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
241
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
242 switch boundary
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
243 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
244 j = 1;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
245 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
246 j = 2;
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
247 otherwise
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
248 error('No such boundary: boundary = %s',boundary);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
249 end
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
250
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
251 switch op
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
252 case 'e'
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
253 switch boundary
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
254 case {'w','W','west','West','s','S','south','South'}
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
255 return_op = obj.e_l{j};
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
256 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
257 return_op = obj.e_r{j};
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
258 end
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
259 case 'd'
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
260 switch boundary
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
261 case {'w','W','west','West','s','S','south','South'}
689
60eb7f46d8d9 Bugfix get_boundary_op in Elastic2d and Heat2d
Martin Almquist <malmquist@stanford.edu>
parents: 688
diff changeset
262 return_op = obj.d1_l{j};
688
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
263 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
689
60eb7f46d8d9 Bugfix get_boundary_op in Elastic2d and Heat2d
Martin Almquist <malmquist@stanford.edu>
parents: 688
diff changeset
264 return_op = obj.d1_r{j};
688
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
265 end
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
266 otherwise
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
267 error(['No such operator: operatr = ' op]);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
268 end
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
269
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
270 end
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
271
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
272 function N = size(obj)
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
273 N = prod(obj.m);
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
274 end
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
275 end
eb2f9233acc3 Add scheme Heat2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
276 end