function [eleidx,mesh,VE]=periodicMesh(num)
% [0,1]^3 
nele=num^3;
nnod=(num+1)^3;
nodenrs=reshape(1:nnod,1+num,1+num,1+num);
edofvec=reshape(nodenrs(1:end-1,1:end-1,1:end-1),nele,1);
edofMat=repmat(edofvec,1,8)+repmat([0 1 num+1 num+2 (num+1)^2+[0 1 num+1 num+2]],nele,1);

nnp=nele;%total number of unique nodes
nnpArray=reshape(1:nnp,num,num,num);
eleidx=nnpArray;%element indices

nnpArray(end+1,:,:)=nnpArray(1,:,:);
nnpArray(:,end+1,:)=nnpArray(:,1,:);
nnpArray(:,:,end+1)=nnpArray(:,:,1);

dofvector=zeros(nnod,1);
dofvector(:)=nnpArray(:);
% dofvector(2:3:end)=3*nnpArray(:)-1;
% dofvector(3:3:end)=3*nnpArray(:);

mesh=dofvector(edofMat);

% coordinates of the 1st node of each element
h=1.0/num;
%V=zeros(num,num,num,3);
VE=zeros(nele,3);
for i=1:num
    for j=1:num
        for k=1:num
            %V(i,j,k,:)=[(i-1)*h,(j-1)*h,(k-1)*h];
            VE(eleidx(i,j,k),:)=[(i-1)*h,(j-1)*h,(k-1)*h];
        end
    end
end
end