Mercurial > repos > public > sbplib

comparison diracDiscrCurve.m @ 1240:1fbd93f24bed feature/dirac_discr

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Factor out Jacobian computation in diracDiscrCurve.m
author Martin Almquist <malmquist@stanford.edu>
date Tue, 19 Nov 2019 16:59:11 -0800
parents 8a456f6e54cc
children 745dc1f1a069
comparison
equal deleted inserted replaced
1238:dea852e85b77 1240:1fbd93f24bed
1 function d = diracDiscrCurve(x_s, g, m_order, s_order, order, opSet) 1 function d = diracDiscrCurve(x_s, g, m_order, s_order, order, opSet)
2 % 2-dimensional delta function for single-block curvilinear grid 2 % 2-dimensional delta function for single-block curvilinear grid
3 % x_s: source point coordinate vector, e.g. [x, y] or [x, y, z]. 3 % x_s: source point coordinate vector, e.g. [x; y] or [x; y; z].
4 % g: single-block grid containing the source 4 % g: single-block grid containing the source
5 % m_order: Number of moment conditions 5 % m_order: Number of moment conditions
6 % s_order: Number of smoothness conditions 6 % s_order: Number of smoothness conditions
7 % order: Order of SBP derivative approximations 7 % order: Order of SBP derivative approximations
8 % opSet: Cell array of function handle to opSet generator 8 % opSet: Cell array of function handle to opSet generator
12 12
13 dim = length(x_s); 13 dim = length(x_s);
14 assert(dim == 2, 'diracDiscrCurve: Only implemented for 2d.'); 14 assert(dim == 2, 'diracDiscrCurve: Only implemented for 2d.');
15 assert(isa(g, 'grid.Curvilinear')); 15 assert(isa(g, 'grid.Curvilinear'));
16 16
17 % Compute Jacobian
18 J = jacobian(g, opSet, order);
19
20 % Find approximate logical coordinates of point source
21 X = g.points();
22 U = g.logic.points();
23 U_interp = scatteredInterpolant(X, U(:,1));
24 V_interp = scatteredInterpolant(X, U(:,2));
25 uS = U_interp(x_s);
26 vS = V_interp(x_s);
27
28 % Get quadrature matrices for moment conditions
29 m = g.size();
30 ops_u = opSet{1}(m(1), {0, 1}, order);
31 ops_v = opSet{2}(m(2), {0, 1}, order);
32 H_u = ops_u.H;
33 H_v = ops_v.H;
34
35 % Get delta function for logical grid and scale by Jacobian
36 d = (1./J).*diracDiscr(g, [uS; vS], m_order, s_order, {H_u, H_v});
37
38 end
39
40 function J = jacobian(g, opSet, order)
17 m = g.size(); 41 m = g.size();
18 m_u = m(1); 42 m_u = m(1);
19 m_v = m(2); 43 m_v = m(2);
20 ops_u = opSet{1}(m_u, {0, 1}, order); 44 ops_u = opSet{1}(m_u, {0, 1}, order);
21 ops_v = opSet{2}(m_v, {0, 1}, order); 45 ops_v = opSet{2}(m_v, {0, 1}, order);
22 I_u = speye(m_u); 46 I_u = speye(m_u);
23 I_v = speye(m_v); 47 I_v = speye(m_v);
24 48
25 D1_u = ops_u.D1; 49 D1_u = ops_u.D1;
26 H_u = ops_u.H;
27
28 D1_v = ops_v.D1; 50 D1_v = ops_v.D1;
29 H_v = ops_v.H;
30 51
31 Du = kr(D1_u,I_v); 52 Du = kr(D1_u,I_v);
32 Dv = kr(I_u,D1_v); 53 Dv = kr(I_u,D1_v);
33 54
34 u = ops_u.x;
35 v = ops_v.x;
36
37 % Compute Jacobian
38 coords = g.points(); 55 coords = g.points();
39 x = coords(:,1); 56 x = coords(:,1);
40 y = coords(:,2); 57 y = coords(:,2);
41 58
42 x_u = Du*x; 59 x_u = Du*x;
43 x_v = Dv*x; 60 x_v = Dv*x;
44 y_u = Du*y; 61 y_u = Du*y;
45 y_v = Dv*y; 62 y_v = Dv*y;
46 63
47 J = x_u.*y_v - x_v.*y_u; 64 J = x_u.*y_v - x_v.*y_u;
48
49 % Find approximate logical coordinates of point source
50 [U, V] = meshgrid(u, v);
51 U_interp = scatteredInterpolant(coords, U(:));
52 V_interp = scatteredInterpolant(coords, V(:));
53 uS = U_interp(x_s);
54 vS = V_interp(x_s);
55
56 d = (1./J).*diracDiscr([uS, vS], {u, v}, m_order, s_order, {H_u, H_v});
57
58 end 65 end