% init 9 kinds of forces
function [center, border, inner, testNum, testF] = forceSamplingInit(nelx,nely,nodenrs)
    
    ix = nelx/2+1; iy = nely/2+1;
    center = (ix-1)*(nely+1)*2+iy*2 + [-1, 0];
    center1 = (ix-1)*(nely+1)*2+(iy+1)*2 + [-1, 0];
    center2 = (ix-1)*(nely+1)*2+(iy-1)*2 + [-1, 0];
    alldofs = 2 * (nely + 1) * (nelx + 1);    
    
    testNum = 3;
    testF = zeros(alldofs, testNum);
    vf = .1;
    
    % four corner points, 1-2-3-4 as usual
    ncorner = [nodenrs(end, [1, end]), nodenrs(1, [end, 1])];
    dcorner = reshape([2 * ncorner - 1; 2 * ncorner], 1, 8); 
    
    % four borders, including corners
    nleft   = nodenrs(1:end, 1)';   dleft   = [2 * nleft - 1, 2 * nleft];
    nbottom = nodenrs(end, 1:end);  dbottom = [2 * nbottom - 1, 2 * nbottom];
    nright  = nodenrs(1:end, end)'; dright  = [2 * nright - 1, 2 * nright];
    nup     = nodenrs(1, 1:end);    dup     = [2 * nup - 1, 2 * nup];
    
    % border points
    nlu = nodenrs(1, 1);   nu = nodenrs(1, (end+1)/2);   nru = nodenrs(1, end);
    nl  = nodenrs((end+1)/2, 1);                         nr  = nodenrs((end+1)/2, end);
    nlb = nodenrs(end, 1); nb = nodenrs(end, (end+1)/2); nrb = nodenrs(end, end);
    
    border = [dcorner, dleft, dbottom, dright, dup];
    
    test = [2*nleft,2*nright];
    inner = setdiff(1:alldofs, [center,center1,center2]);
    
    
%     testF(2 * [nlu, nlb] - 1, j) = -vf/4;
%     testF(2 * [nru, nrb] - 1, j) = vf/4;
%     testF(2 * [nl] - 1, j) = -vf;
%     testF(2 * [nr] - 1, j) = vf;
    % horizontal tensile - 10
    testF(2 * nleft - 1, 1) = -vf;
    testF(2 * nright - 1, 1) = vf;
%     for j = 1:10
%         testF(2 * nleft - 1, j) = -vf*(j-1)/9;
%         testF(2 * nright - 1, j) = vf*(1-(j-1)/9);
%     end
        
%     testF(2 * [nlu, nru], j) = vf/4;
%     testF(2 * [nlb, nrb], j) = -vf/4;
%     testF(2 * [nu], j) = vf;
%     testF(2 * [nb], j) = -vf;
    % longitudinal tensile - 10
    testF(2 * nup, 2) = vf;
    testF(2 * nbottom, 2) = -vf;
%     for j = 11:20
%         testF(2 * nup, j) = vf*(1-(j-11)/9);
%         testF(2 * nbottom, j) = -vf*(j-11)/9;
%     end
        
%     testF(2 * [nu, nru,nr]-1, j) = vf;
%     testF(2 * [nu, nru,nr], j) = vf;
%     testF(2 * [nl, nlb, nb]-1, j) = -vf;
%     testF(2 * [nl, nlb, nb], j) = -vf;    
    % shear - 20
    testF([2 * nup-1, 2 * nright - 1], 3) = vf/2;
    testF([2 * nup, 2 * nright], 3) = vf/2;
    testF([2 * nleft-1, 2 * nbottom - 1], 3) = -vf/2;
    testF([2 * nleft, 2 * nbottom], 3) = -vf/2; 
    testF([2*[nlu,nrb]-1, 2*[nlu,nrb]], 3) = 0;
%     for j = (21:30)
%         testF([2 * nup-1, 2 * nright - 1], j) = (1-(j-21)/9)*vf/2;
%         testF([2 * nup, 2 * nright], j) = (1-(j-21)/9)*vf/2;
%         testF([2 * nleft-1, 2 * nbottom - 1], j) = -vf/2*(j-21)/9;
%         testF([2 * nleft, 2 * nbottom], j) = -vf/2*(j-21)/9; 
%         testF([2*[nlu,nrb]-1, 2*[nlu,nrb]], j) = 0;
%     end
    
    j = 4; % bending - 20
    
    j = 5; % twisting - 20
    
    j = 6; % pinching - 20
end