Mercurial > repos > public > sbplib

annotate +multiblock/LaplaceSquared.m @ 1301:8978521b0f06 default

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix incorrect package name.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 08 Jul 2020 19:11:04 +0200
parents 2aaced07d1e5
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
1070
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
1 classdef LaplaceSquared < scheme.Scheme
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
2 properties
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
3 grid
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
4 order
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
5 laplaceDiffOp
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
6
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
7 D
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
8 H
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
9 Hi
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
10
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
11 a,b
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
12 end
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
13
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
14 methods
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
15 % Discretisation of a*nabla*b*nabla
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
16 function obj = LaplaceSquared(g, order, a, b, opGen)
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
17 default_arg('order', 4);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
18 default_arg('a', 1);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
19 default_arg('b', 1);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
20 default_arg('opGen', @sbp.D4Variable);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
21
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
22 if isscalar(a)
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
23 a = grid.evalOn(g, a);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
24 end
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
25
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
26 if isscalar(b)
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
27 b = grid.evalOn(g, b);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
28 end
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
29
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
30 obj.grid = g;
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
31 obj.order = order;
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
32 obj.a = a;
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
33 obj.b = b;
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
34
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
35 obj.laplaceDiffOp = multiblock.Laplace(g, order, 1, 1, opGen);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
36
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
37 obj.H = obj.laplaceDiffOp.H;
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
38 obj.Hi = spdiag(1./diag(obj.H));
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
39
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
40 A = spdiag(a);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
41 B = spdiag(b);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
42
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
43 D_laplace = obj.laplaceDiffOp.D;
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
44 obj.D = A*D_laplace*B*D_laplace;
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
45 end
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
46
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
47 function s = size(obj)
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
48 s = size(obj.laplaceDiffOp);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
49 end
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
50
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
51 function op = getBoundaryOperator(obj, opName, boundary)
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
52 switch opName
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
53 case 'e'
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
54 op = getBoundaryOperator(obj.laplaceDiffOp, 'e', boundary);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
55 case 'd1'
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
56 op = getBoundaryOperator(obj.laplaceDiffOp, 'd', boundary);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
57 case 'd2'
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
58 e = getBoundaryOperator(obj.laplaceDiffOp, 'e', boundary);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
59 op = (e'*obj.laplaceDiffOp.D)';
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
60 case 'd3'
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
61 d1 = getBoundaryOperator(obj.laplaceDiffOp, 'd', boundary);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
62 op = (d1'*spdiag(obj.b)*obj.laplaceDiffOp.D)';
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
63 end
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
64 end
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
65
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
66 function op = getBoundaryQuadrature(obj, boundary)
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
67 op = getBoundaryQuadrature(obj.laplaceDiffOp, boundary);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
68 end
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
69
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
70 function [closure, penalty] = boundary_condition(obj,boundary,type) % TODO: Change name to boundaryCondition
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
71 switch type
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
72 case 'e'
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
73 error('Bc of type ''e'' not implemented')
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
74 case 'd1'
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
75 error('Bc of type ''d1'' not implemented')
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
76 case 'd2'
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
77 e = obj.getBoundaryOperator('e', boundary);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
78 d1 = obj.getBoundaryOperator('d1', boundary);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
79 d2 = obj.getBoundaryOperator('d2', boundary);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
80 H_b = obj.getBoundaryQuadrature(boundary);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
81
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
82 A = spdiag(obj.a);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
83 B_b = spdiag(e'*obj.b);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
84
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
85 tau = obj.Hi*A*d1*B_b*H_b;
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
86 closure = tau*d2';
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
87 penalty = -tau;
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
88 case 'd3'
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
89 e = obj.getBoundaryOperator('e', boundary);
1073
95113a592421 Fix bug in d3 boundary condition
Jonatan Werpers <jonatan@werpers.com>
parents: 1070
diff changeset
90 d3 = obj.getBoundaryOperator('d3', boundary);
1070
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
91 H_b = obj.getBoundaryQuadrature(boundary);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
92
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
93 A = spdiag(obj.a);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
94
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
95 tau = -obj.Hi*A*e*H_b;
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
96 closure = tau*d3';
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
97 penalty = -tau;
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
98 end
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
99 end
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
100
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
101 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
102 error('Not implemented')
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
103 end
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
104 end
1074
2aaced07d1e5 Formatting fix
Jonatan Werpers <jonatan@werpers.com>
parents: 1073
diff changeset
105 end