2 prepNodesAndWeights(tn,dt,beta,gl,s,src_dmn,snIdxNow,P,M,W,order,...
3 timeLevels,typNumOfNeighbors,winData)
4% PREPNODESANDWEIGHTS to prepare nodes and weights for the numerical
5% evaluation of integrals
6%
7% INPUTS:
8% tn,dt: time grid (array) with time step dt
9% phi_wts: array of weights to evaluate the blending function using
10% Chebyshev polynomial
11% beta: an array containing spring constants
12% gl: GL nodes gl.x0 and weights gl.w0 over the interval [-1,1];
13% s: a vector holding the locations of the spring sources
14% snIdxNow: index of density grid function at the current step. For
15% example, the density functions are arranged in the following order
16% -tgh, ..., -1,0,1, where tgh corresponds to temporal ghost points,
17% snIdxNow corresponds to the index of 0.
18% P: number of GL nodes
19% M: number of sources
20% W: width of the window function (i.e. support of the window function up
21% to tol)
22% chebApproxInfo: information needed for the Chebyshev polynomial
23% approximation of window and blending functions
24% timeLevels: the number of time levels stored for the density functions;
25% this is equal to W + order to account for centered barycentric
26% interpolation at the left end of the time axis
27% typNumOfNeighbors: estimated average typical number of neighbors for
28% spring sources
29%
30% OUTPUTS:
31% implicitMat: implicit sparse matrix to solve for the density values at
32% the new time
33% A_sp: a sparse matrix to evaluate integrals numerically using a
34% combination of barycentric interpolation and Gauss-Legendre quadrature
36% domain of sources
37as = src_dmn(1); bs = src_dmn(2);
39% support of the window function
40delta = W*dt;
42% prepare the implicit matrix
43implicitMat = speye(M);
45% get the LMM weights for the treatment of selfLocal1
46selfLocal1_wts = get_lmmwts(tn(1:order), order, order,dt,dt,delta,winData,0);
47coef = -(1 + (selfLocal1_wts(end).*(beta./2))); % coefficient of the sigma at (n+1)
48implicitMat = coef.*implicitMat;
50% GL nodes and weights for selfLocal2 and blending function on [dt,delta]
51[eta,~] = glwt(dt,delta,gl);
52phidt_to_delta = generalwindow((eta./delta),winData);
54% GL nodes and weights for the selfLocal2 implementation and
55% interpolation weights
56[selfLocal2_glnodes,selfLocal2_glwts] = glwt(0,(delta-dt),gl);
57selfLocal2_interpWeights_temp = zeros(P,order);
58interpNodesIdx_temp = zeros(P,order);
59for kk = 1:P % for each GL node calculate interpolated value at GL node
60 [interpNodes,interpNodesIdx_temp(kk,:)] = get_InterpNodes(selfLocal2_glnodes(kk), dt, order,(delta-dt));
61 selfLocal2_interpWeights_temp(kk,:) = barycentricInterp(interpNodes, order, selfLocal2_glnodes(kk));
62end
63% multiply the weights by GL weights and window factor
64selfLocal2_interpWeights_temp = selfLocal2_interpWeights_temp.*(1 - phidt_to_delta(end:-1:1)).*selfLocal2_glwts;
66% shift the nodes to fit the position of each node
67interpNodesIdx_temp = interpNodesIdx_temp+ snIdxNow + 1 - W;
69% merge the weights of similar interpolation nodes
70[selfLocal2_interpWeights] ...
71 = mergeWeights(interpNodesIdx_temp,selfLocal2_interpWeights_temp,timeLevels);
73% create bins to avoid checks over the full set of sources
74% each box is delta big, so check sources in the current box and other
75% surrounding boxes
76nboxes = ceil((bs - as)/delta); % number of boxes
77snew = (s - as)./(bs - as); % transform s so it is on [0,1] for assign function
78[ioffst, ibox, isradr,icnt] = assign(nboxes, snew, M); % assign sources to boxes
80% Fill in a sparse matrix with weights to evaluate within the time loop
81% The product of the sparse matrix with the density vector evaluates
82% numerical integrals needed in the solver
83nz = M*timeLevels*typNumOfNeighbors;
84i_sp = ones(1,nz);
85j_sp = ones(1,nz);
86a_sp = zeros(1,nz);
88cnt = 1;
89dtInteractions = zeros(M,1); % set up counts for delta interactions for testing
90deltaInteractions = zeros(M,1); % set up counts for dt interactions for testing
91for j = 1:M
92 adr = ibox(j);
93 sourcesNearby_left = []; sourcesNearby_right = [];
94 sourcesNearby_center = isradr(ioffst(adr):(ioffst(adr) + icnt(adr) - 1));
95 if(adr>1)
96 sourcesNearby_left = [isradr(ioffst(adr-1):(ioffst(adr-1) + icnt(adr-1) - 1))];
97 end
98 if(adr<nboxes)
99 sourcesNearby_right = [isradr(ioffst(adr+1):(ioffst(adr+1) + icnt(adr+1) - 1))];
100 end
101 sourcesNearby = [sourcesNearby_center; sourcesNearby_left; sourcesNearby_right];
102 for l = sourcesNearby'
103 % prepare the indices of the sparse matrix
104 lInd = ((l-1)*timeLevels + 1):((l-1)*timeLevels + timeLevels);
105 Idx = ((cnt - 1)*timeLevels+1):((cnt - 1)*timeLevels+timeLevels);
107 if(l~=j)
108 % local evaluation from other sources
109 L = abs(s(j) - s(l));
110 if(L<delta && L>dt)
111 % get GL nodes and weights for the integral involving other
112 % neighboring sources, but does not include the density at
113 % the new step
114 [otherLocal2_glnodes,otherLocal2_glwts] = glwt(0,(delta - L),gl);
115 otherLocal2_interpWeights_temp = zeros(P,order);
116 otherInterpNodesIdx_temp = zeros(P,order);
117 for kk = 1:P % for each GL node calculate interpolated value at GL node
118 [otherInterpNodes,otherInterpNodesIdx_temp(kk,:)] = get_InterpNodes(otherLocal2_glnodes(kk), dt, order,(delta - L));
119 otherLocal2_interpWeights_temp(kk,:) = barycentricInterp(otherInterpNodes, order, otherLocal2_glnodes(kk));
120 end
121 % prepare the values of phi to evaluate the local
122 % contribution from each of the other sources
123 [eta,~] = glwt(L,delta,gl);
124 phiL_to_delta = generalwindow((eta./delta),winData);
126 % multiply the interpolation weights with the blending
127 % function
128 otherLocal2_interpWeights_temp = otherLocal2_interpWeights_temp.*(1 - phiL_to_delta(end:-1:1)).*otherLocal2_glwts;
130 % shift the indices;
131 otherInterpNodesIdx_temp = otherInterpNodesIdx_temp + snIdxNow + 1 - W;
133 [otherLocal2_interpWeights] ...
134 = mergeWeights(otherInterpNodesIdx_temp,otherLocal2_interpWeights_temp,timeLevels);
136 i_sp(Idx) = j;
137 j_sp(Idx) = lInd;
138 a_sp(Idx) = otherLocal2_interpWeights;
140 cnt = cnt + 1 ;
141 deltaInteractions(j) = deltaInteractions(j) + 1;
143 elseif(L<dt)
145 otherLocal1_wts = get_lmmwts(tn(1:order), order, order,(dt - L),dt,delta,winData,0);
146 coef_other = (-beta(j)/2)*otherLocal1_wts(end); % coefficient of the sigma at (n+1)
147 implicitMat(j,l) = coef_other;
149 i_sp(Idx) = j;
150 j_sp(Idx) = lInd;
151 a_sp(Idx) = selfLocal2_interpWeights;
153 % Adding the local 1 contribution
154 endTerm = Idx(end)-1; % Index refering tothe current step
155 a_sp((endTerm+1-(order-1)):endTerm) = a_sp((endTerm+1-(order-1)):endTerm) + otherLocal1_wts((end-(order-1)):end-1);
156 cnt = cnt + 1 ;
158 dtInteractions(j) = dtInteractions(j) + 1;
159 end
160 elseif(j == l)
161 % Adding the self local 2 contributions
162 i_sp(Idx) = j;
163 j_sp(Idx) = lInd;
164 a_sp(Idx) = selfLocal2_interpWeights;
166 % Adding the self local 1 contributions
167 endTerm = Idx(end)-1; % Index refering to the current step
168 a_sp((endTerm+1-(order-1)):endTerm) = a_sp((endTerm+1-(order-1)):endTerm) + selfLocal1_wts((end-(order-1)):end-1);
169 cnt = cnt + 1 ;
170 end
172 end % end of l loop
173end % end of j loop
175% decompose the matrix to speed up solution in the time-stepping loop
176implicitMat = decomposition(implicitMat);
178A_sp = sparse(i_sp,j_sp,a_sp,M,(M*timeLevels));
180clear i_sp j_sp a_sp;
182end % end of function