8748c66added codenalhassanieh 1function ul = get_uLocal(x,t,snIdxNow,sn,interpNodesIdx_sol,interpAndGLWeights_sol,interpShift_sol,s,M,delta)
2% GET_ULOCAL computes the local part of the solution
3%
4% ul = get_uLocal(x,t,snIdxNow,sn,interpNodesIdx_sol,...
5% interpAndGLWeights_sol,interpShift_sol,s,M,delta)
6% returns the local part of the solution at each value of the vector 'x' at
7% time 't'. The function takes the unifrom time-grid 'tn', 'delta' ([0,delta]
8% is the compact support of the window), 'M' source locations in the vector
9% s = [s(1),...,s(M)], the computed density values sn.
10% snIdxNow: index of the density at the current time step.
11% interpNodesIdx_sol: interpolation nodes used to evaluate the density at
12% GL nodes
13% interpAndGLWeights_sol: interpolation/integration weights needed to
14% evaluate the local integral
15% interpShift_sol: initial shift of interpolation nodes (to undo the shift)
17% Pick interpolation points between [startTime,endTime], where startTime
18% and endTime lie on the time-grid
19Nx = length(x);
20 ul = zeros(Nx,1);
21 for i = 1:Nx % for each x value
22 for j = 1:M % for each source
23 L = abs(x(i) - s(j));
24 if(L<delta && (t - L)>0)
25 snj = sn(:,j);
26 % perform the local integration using GL and barycentric
27 % interpolation
28 shift = (snIdxNow - 1) + interpShift_sol;
30 shiftedInterpNodesIdx = interpNodesIdx_sol{i,j} + shift;
31 snIdx = shiftedInterpNodesIdx + 1; %snIdx(snIdx<=0) = 1;
32 snvec = sum(interpAndGLWeights_sol{i,j}.*snj(snIdx),2);
34 I = sum(snvec);
35 ul(i) = ul(i) + I;
36 end
37 end
38 ul(i) = 0.5*ul(i);
39 end
41end
43% Note replaced for loops by matvecs. Original:
44% snvec = zeros(P,1);
45% for k = 1:P % for each GL node calculate interpolated value at GL node
46% shiftedInterpNodesIdx = interpNodesIdx_sol{i,j}(k,:) + shift;
47% snIdx = shiftedInterpNodesIdx + 1; snIdx(snIdx<=0) = 1;
48% snvec(k) = interpAndGLWeights_sol{i,j}(k,:)*(sn{j}(snIdx));
49% end
50% I = sum(snvec);