Mercurial > repos > public > sbplib

comparison diracDiscr.m @ 1241:d1b201fe328e feature/dirac_discr

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Turn row vectors into column vectors in diracDiscr comments
author Martin Almquist <malmquist@stanford.edu>
date Tue, 19 Nov 2019 17:01:38 -0800
parents dea852e85b77
children ff613067dec6
comparison
equal deleted inserted replaced
1240:1fbd93f24bed 1241:d1b201fe328e
1 1
2 function d = diracDiscr(g, x_s, m_order, s_order, H) 2 function d = diracDiscr(g, x_s, m_order, s_order, H)
3 % n-dimensional delta function 3 % n-dimensional delta function
4 % g: cartesian grid 4 % g: cartesian grid
5 % x_s: source point coordinate vector, e.g. [x, y] or [x, y, z]. 5 % x_s: source point coordinate vector, e.g. [x; y] or [x; y; z].
6 % m_order: Number of moment conditions 6 % m_order: Number of moment conditions
7 % s_order: Number of smoothness conditions 7 % s_order: Number of smoothness conditions
8 % H: cell array of 1D norm matrices 8 % H: cell array of 1D norm matrices
9 assertType(g, 'grid.Cartesian'); 9 assertType(g, 'grid.Cartesian');
10 dim = g.d; 10 dim = g.d;
28 28
29 end 29 end
30 30
31 31
32 % Helper function for 1D delta functions 32 % Helper function for 1D delta functions
33 function ret = diracDiscr1D(x_s , x , m_order, s_order, H) 33 function ret = diracDiscr1D(x_s, x, m_order, s_order, H)
34 34
35 m = length(x); 35 m = length(x);
36 36
37 % Return zeros if x0 is outside grid 37 % Return zeros if x0 is outside grid
38 if x_s < x(1) || x_s > x(end) 38 if x_s < x(1) || x_s > x(end)
49 49
50 index = sourceIndices(x_s, x, tot_order, h); 50 index = sourceIndices(x_s, x, tot_order, h);
51 51
52 polynomial = (x(index)-x(index(1)))/(x(index(end))-x(index(1))); 52 polynomial = (x(index)-x(index(1)))/(x(index(end))-x(index(1)));
53 x_0 = (x_s-x(index(1)))/(x(index(end))-x(index(1))); 53 x_0 = (x_s-x(index(1)))/(x(index(end))-x(index(1)));
54 54
55 quadrature = diag(H); 55 quadrature = diag(H);
56 quadrature_weights = quadrature(index)/h; 56 quadrature_weights = quadrature(index)/h;
57 57
58 h_polynomial = polynomial(2)-polynomial(1); 58 h_polynomial = polynomial(2)-polynomial(1);
59 b = zeros(tot_order,1); 59 b = zeros(tot_order,1);