Mercurial > repos > public > sbplib

annotate diracDiscrCurve.m @ 1336:0666629aa183 feature/D2_boundary_opt

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add methods for creating grids with different grid point distributions for each coordinate direction, and also supports constructing periodic grids
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 13 May 2022 13:26:16 +0200
parents 25efceb0c392
children 663eb90a4559
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
1246
25efceb0c392 Apply some nitpicking
Jonatan Werpers <jonatan@werpers.com>
parents: 1244
diff changeset
1 % 2-dimensional delta function for single-block curvilinear grid
25efceb0c392 Apply some nitpicking
Jonatan Werpers <jonatan@werpers.com>
parents: 1244
diff changeset
2 % x_s: source point coordinate vector, e.g. [x; y] or [x; y; z].
25efceb0c392 Apply some nitpicking
Jonatan Werpers <jonatan@werpers.com>
parents: 1244
diff changeset
3 % g: single-block grid containing the source
25efceb0c392 Apply some nitpicking
Jonatan Werpers <jonatan@werpers.com>
parents: 1244
diff changeset
4 % m_order: Number of moment conditions
25efceb0c392 Apply some nitpicking
Jonatan Werpers <jonatan@werpers.com>
parents: 1244
diff changeset
5 % s_order: Number of smoothness conditions
25efceb0c392 Apply some nitpicking
Jonatan Werpers <jonatan@werpers.com>
parents: 1244
diff changeset
6 % order: Order of SBP derivative approximations
25efceb0c392 Apply some nitpicking
Jonatan Werpers <jonatan@werpers.com>
parents: 1244
diff changeset
7 % opSet: Cell array of function handle to opSet generator
1230
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
8 function d = diracDiscrCurve(x_s, g, m_order, s_order, order, opSet)
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
9 default_arg('order', m_order);
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
10 default_arg('opSet', {@sbp.D2Variable, @sbp.D2Variable});
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
11
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
12 dim = length(x_s);
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
13 assert(dim == 2, 'diracDiscrCurve: Only implemented for 2d.');
1244
745dc1f1a069 Use assertType instead of assert(isa(...))
Jonatan Werpers <jonatan@werpers.com>
parents: 1240
diff changeset
14 assertType(g, 'grid.Curvilinear');
1230
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
15
1240
1fbd93f24bed Factor out Jacobian computation in diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents: 1230
diff changeset
16 % Compute Jacobian
1fbd93f24bed Factor out Jacobian computation in diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents: 1230
diff changeset
17 J = jacobian(g, opSet, order);
1fbd93f24bed Factor out Jacobian computation in diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents: 1230
diff changeset
18
1fbd93f24bed Factor out Jacobian computation in diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents: 1230
diff changeset
19 % Find approximate logical coordinates of point source
1fbd93f24bed Factor out Jacobian computation in diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents: 1230
diff changeset
20 X = g.points();
1fbd93f24bed Factor out Jacobian computation in diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents: 1230
diff changeset
21 U = g.logic.points();
1fbd93f24bed Factor out Jacobian computation in diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents: 1230
diff changeset
22 U_interp = scatteredInterpolant(X, U(:,1));
1fbd93f24bed Factor out Jacobian computation in diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents: 1230
diff changeset
23 V_interp = scatteredInterpolant(X, U(:,2));
1fbd93f24bed Factor out Jacobian computation in diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents: 1230
diff changeset
24 uS = U_interp(x_s);
1fbd93f24bed Factor out Jacobian computation in diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents: 1230
diff changeset
25 vS = V_interp(x_s);
1fbd93f24bed Factor out Jacobian computation in diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents: 1230
diff changeset
26
1fbd93f24bed Factor out Jacobian computation in diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents: 1230
diff changeset
27 % Get quadrature matrices for moment conditions
1fbd93f24bed Factor out Jacobian computation in diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents: 1230
diff changeset
28 m = g.size();
1fbd93f24bed Factor out Jacobian computation in diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents: 1230
diff changeset
29 ops_u = opSet{1}(m(1), {0, 1}, order);
1fbd93f24bed Factor out Jacobian computation in diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents: 1230
diff changeset
30 ops_v = opSet{2}(m(2), {0, 1}, order);
1fbd93f24bed Factor out Jacobian computation in diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents: 1230
diff changeset
31 H_u = ops_u.H;
1fbd93f24bed Factor out Jacobian computation in diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents: 1230
diff changeset
32 H_v = ops_v.H;
1fbd93f24bed Factor out Jacobian computation in diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents: 1230
diff changeset
33
1fbd93f24bed Factor out Jacobian computation in diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents: 1230
diff changeset
34 % Get delta function for logical grid and scale by Jacobian
1fbd93f24bed Factor out Jacobian computation in diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents: 1230
diff changeset
35 d = (1./J).*diracDiscr(g, [uS; vS], m_order, s_order, {H_u, H_v});
1fbd93f24bed Factor out Jacobian computation in diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents: 1230
diff changeset
36 end
1fbd93f24bed Factor out Jacobian computation in diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents: 1230
diff changeset
37
1fbd93f24bed Factor out Jacobian computation in diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents: 1230
diff changeset
38 function J = jacobian(g, opSet, order)
1230
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
39 m = g.size();
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
40 m_u = m(1);
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
41 m_v = m(2);
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
42 ops_u = opSet{1}(m_u, {0, 1}, order);
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
43 ops_v = opSet{2}(m_v, {0, 1}, order);
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
44 I_u = speye(m_u);
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
45 I_v = speye(m_v);
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
46
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
47 D1_u = ops_u.D1;
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
48 D1_v = ops_v.D1;
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
49
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
50 Du = kr(D1_u,I_v);
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
51 Dv = kr(I_u,D1_v);
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
52
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
53 coords = g.points();
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
54 x = coords(:,1);
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
55 y = coords(:,2);
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
56
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
57 x_u = Du*x;
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
58 x_v = Dv*x;
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
59 y_u = Du*y;
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
60 y_v = Dv*y;
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
61
8a456f6e54cc Add diracDiscrCurve.m
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
62 J = x_u.*y_v - x_v.*y_u;
1246
25efceb0c392 Apply some nitpicking
Jonatan Werpers <jonatan@werpers.com>
parents: 1244
diff changeset
63 end