1function [interpNodesIdx_sol,interpAndGLWeights_sol,interpShift_sol] = prepForLocalSolEval(x,dt,gl,s,M,P,W,order,winData)
2% PREPFORLOCALSOLEVAL precomputes parameters for local solution evaluation
3%
4% [interpNodesIdx_sol,interpAndGLWeights_sol,interpShift_sol] = ...
5% prepForLocalSolEval(x,dt,delta,gl,s,M,phi_wts,P,W,order,chebApproxInfo)
6% prepares interpolation nodes and weights needed for the evaluation
7% of the local part of the solution. The function takes 'x' the
8% spatial grid, 'dt' the time-step, 'delta = W*dt' the width of the window
9% function 'gl' GL nodes and weights on [-1,1], 's' a vector
10% carrying the location of the sources, 'M' the number of sources, the
11% blending function cheb weights 'phi_wts' and the number of GL nodes 'P'.
12%
13% NOTE:'chebApproxInfo': a struct with fields domain [ax,bx] over which
14% cheb weights are computed; ninters, number of subintervals, and nord,
15% order of polynomial interpolation or number of cheb nodes
16%
17% OUTPUT:
18% interpNodesIdx_sol: indices of interpolation nodes
19% interpAndGLWeights_sol: weights needed to perform the integration of
20% the local part of the soluiton
21% intepShift_sol: how many grid points are the interpolation points
22% shifted. We undo the shift when we use these nodes in the
23% time-stepping loop.
25delta = W*dt;
27Nx = length(x); % number of points in the spatial grid
29% allocate space for variables
30interpNodesIdx_sol = cell(Nx,M); interpNodesIdx_sol(:) = {zeros(P,order)};
31interpAndGLWeights_sol = cell(Nx,M); interpAndGLWeights_sol(:) = {zeros(P,order)};
33interpShift_sol = -(W-1); % shift index of the interpolation nodes (nodes shifted forward by delta)
35for i = 1:Nx % for each x vlue
36 for j = 1:M % for each source
37 L = abs(x(i) - s(j));
38 if(L<delta)
40 % prepare the GL nodes eta needed for window factor
41 [eta,~] = glwt(L,delta,gl);
43 % Get the blending function at eta
44 phiL_to_delta_sol = generalwindow((eta./delta),winData);
46 % actual interval of integration
47 % t1 = t - delta;
48 % t2 = t - L;
50 % interval of integration (will be shifted later)
51 t1 = 0;
52 t2 = delta - L;
54 % Get the GL nodes and weights in the interval of integration
55 [glnodes,glwts] = glwt(t1,t2,gl);
57 % update the glwts by the window factor
58 glwts = glwts.*(1 - phiL_to_delta_sol(end:-1:1));
60 % perform Barycentric interpolation at the GL nodes and get
61 % weights
62 interpWeights_sol = zeros(P,order);
63 endTime = t2;
64 for k = 1:P
65 [interpNodes,interpNodesIdx_sol{i,j}(k,:)] = get_InterpNodes(glnodes(k), dt, order,endTime);
66 interpWeights_sol(k,:) = barycentricInterp(interpNodes, order, glnodes(k));
67 end
69 % Multiply the interpolation weights by the GL weights
70 interpAndGLWeights_sol{i,j} = interpWeights_sol.*glwts;
71 end
72 end
73end
75end