Mercurial > repos > public > sbplib

annotate +scheme/Gradient.m @ 1339:bcdb14b05d03 feature/D2_boundary_opt

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix comment
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 22 Jul 2022 16:31:18 +0200
parents dd25184e71ec
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
1121
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
1 classdef Gradient < scheme.Scheme
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
2
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
3 % Approximates the divergence
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
4 % Interface and boundary condition methods are just dummies
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
5
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
6 properties
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
7 m % Number of points in each direction, possibly a vector
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
8 h % Grid spacing
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
9
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
10 grid
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
11 dim
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
12
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
13 order % Order of accuracy for the approximation
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
14
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
15 D
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
16 D1
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
17 H
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
18 end
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
19
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
20 methods
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
21
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
22 function obj = Gradient(g, order, opSet)
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
23 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
24
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
25 dim = 2;
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
26
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
27 m = g.size();
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
28 m_tot = g.N();
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
29
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
30 h = g.scaling();
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
31 lim = g.lim;
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
32 if isempty(lim)
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
33 x = g.x;
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
34 lim = cell(length(x),1);
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
35 for i = 1:length(x)
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
36 lim{i} = {min(x{i}), max(x{i})};
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
37 end
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
38 end
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
39
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
40 % 1D operators
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
41 ops = cell(dim,1);
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
42 for i = 1:dim
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
43 ops{i} = opSet{i}(m(i), lim{i}, order);
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
44 end
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
45
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
46 I = cell(dim,1);
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
47 D1 = cell(dim,1);
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
48
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
49 for i = 1:dim
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
50 I{i} = speye(m(i));
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
51 D1{i} = ops{i}.D1;
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
52 end
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
53
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
54 %====== Assemble full operators ========
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
55
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
56 % D1
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
57 obj.D1{1} = kron(D1{1},I{2});
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
58 obj.D1{2} = kron(I{1},D1{2});
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
59
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
60 I_dim = speye(dim, dim);
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
61
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
62 % E{i}^T picks out component i.
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
63 E = cell(dim,1);
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
64 I = speye(m_tot,m_tot);
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
65 for i = 1:dim
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
66 e = sparse(dim,1);
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
67 e(i) = 1;
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
68 E{i} = kron(I,e);
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
69 end
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
70
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
71 Grad = sparse(dim*m_tot, m_tot);
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
72 for i = 1:dim
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
73 Grad = Grad + E{i}*obj.D1{i};
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
74 end
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
75 obj.D = Grad;
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
76 obj.H = [];
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
77 %=========================================%'
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
78
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
79 end
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
80
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
81
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
82 % Closure functions return the operators applied to the own domain to close the boundary
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
83 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
84 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
85 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
86 % on the first component. Can also be e.g.
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
87 % {'normal', 'd'} or {'tangential', 't'} for conditions on
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
88 % tangential/normal component.
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
89 % data is a function returning the data that should be applied at the boundary.
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
90 % neighbour_scheme is an instance of Scheme that should be interfaced to.
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
91 % neighbour_boundary is a string specifying which boundary to interface to.
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
92 function [closure, penalty] = boundary_condition(obj, boundary, bc)
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
93 error('Not implemented')
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
94 end
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
95
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
96 % type Struct that specifies the interface coupling.
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
97 % Fields:
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
98 % -- tuning: penalty strength, defaults to 1.2
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
99 % -- interpolation: type of interpolation, default 'none'
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
100 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
101
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
102 [m, n] = size(obj.D);
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
103 closure = sparse(m, n);
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
104
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
105 [m, n] = size(neighbour_scheme.D);
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
106 penalty = sparse(m, n);
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
107 end
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
108
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
109 function N = size(obj)
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
110 N = obj.dim*prod(obj.m);
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
111 end
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
112 end
dd25184e71ec Add diffOp for gradient
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
113 end