|
1230
|
1 function d = diracDiscrCurve(x_s, g, m_order, s_order, order, opSet)
|
|
|
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].
|
|
|
4 % g: single-block grid containing the source
|
|
|
5 % m_order: Number of moment conditions
|
|
|
6 % s_order: Number of smoothness conditions
|
|
|
7 % order: Order of SBP derivative approximations
|
|
|
8 % opSet: Cell array of function handle to opSet generator
|
|
|
9
|
|
|
10 default_arg('order', m_order);
|
|
|
11 default_arg('opSet', {@sbp.D2Variable, @sbp.D2Variable});
|
|
|
12
|
|
|
13 dim = length(x_s);
|
|
|
14 assert(dim == 2, 'diracDiscrCurve: Only implemented for 2d.');
|
|
|
15 assert(isa(g, 'grid.Curvilinear'));
|
|
|
16
|
|
|
17 m = g.size();
|
|
|
18 m_u = m(1);
|
|
|
19 m_v = m(2);
|
|
|
20 ops_u = opSet{1}(m_u, {0, 1}, order);
|
|
|
21 ops_v = opSet{2}(m_v, {0, 1}, order);
|
|
|
22 I_u = speye(m_u);
|
|
|
23 I_v = speye(m_v);
|
|
|
24
|
|
|
25 D1_u = ops_u.D1;
|
|
|
26 H_u = ops_u.H;
|
|
|
27
|
|
|
28 D1_v = ops_v.D1;
|
|
|
29 H_v = ops_v.H;
|
|
|
30
|
|
|
31 Du = kr(D1_u,I_v);
|
|
|
32 Dv = kr(I_u,D1_v);
|
|
|
33
|
|
|
34 u = ops_u.x;
|
|
|
35 v = ops_v.x;
|
|
|
36
|
|
|
37 % Compute Jacobian
|
|
|
38 coords = g.points();
|
|
|
39 x = coords(:,1);
|
|
|
40 y = coords(:,2);
|
|
|
41
|
|
|
42 x_u = Du*x;
|
|
|
43 x_v = Dv*x;
|
|
|
44 y_u = Du*y;
|
|
|
45 y_v = Dv*y;
|
|
|
46
|
|
|
47 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 |
