Faceids for Modflow example

Faceids_4_Modflow_example.m

faceVel=nan(15000,22);
faceVel(1,:) = 0;
cnt_fc = 1;
msh_face_ids = zeros(size(msh,1),4);
cnt = 1;
for ii = 1:size(msh_grid,1)
    for jj = 1:size(msh_grid,2)
        nd = msh_grid(ii,jj).ND;
        % Face 1 is the ii-1,jj FRONT
        if ii == 1
            msh_face_ids(cnt,1) = -1;
        else
            idx = find(faceVel(1:cnt_fc,1) == nd(2) & faceVel(1:cnt_fc,2) == nd(1));
            if isempty(idx)
                idx = find(faceVel(1:cnt_fc,1) == nd(1) & faceVel(1:cnt_fc,2) == nd(2));
            end
            if isempty(idx)
                stop_here = true;
            else
                msh_face_ids(cnt,1) = -idx;
            end
        end
        
        % Face 2 is the ii,jj Right
        if jj ~= 120
            cnt_fc = cnt_fc + 1;
            faceVel(cnt_fc,:) = [nd(2:3) msh_grid(ii,jj).RIGHT'];
            msh_face_ids(cnt,2) = cnt_fc;
        end
        
        
        % Face 3 is the ii,jj FRONT
        if ii == 60
            msh_face_ids(cnt,3) = 1;
        else
            cnt_fc = cnt_fc + 1;
            faceVel(cnt_fc,:) = [nd(3:4) msh_grid(ii,jj).FRONT'];
            msh_face_ids(cnt,3) = cnt_fc;
        end
        
        % Face 4 is the is the ii,jj-1 RIGHT but we will use the ii,jj
        % RIGHT for the first column
        if jj == 1 
            msh_face_ids(cnt,4) = -msh_face_ids(cnt,2);
        else
            idx = find(faceVel(1:cnt_fc,1) == nd(1) & faceVel(1:cnt_fc,2) == nd(4));
            if isempty(idx)
                idx = find(faceVel(1:cnt_fc,1) == nd(4) & faceVel(1:cnt_fc,2) == nd(1));
            end
            if isempty(idx)
                stop_here = true;
            else
                msh_face_ids(cnt,4) = -idx;
                if jj == 120
                    msh_face_ids(cnt,2) = idx;
                end
            end
        end
        
        cnt = cnt + 1;
    end
end
faceVel(cnt_fc+1:end,:) = [];