function CreateGradBasisFunction3D(iscale,nmod)
%sbasis,frestart
load(fullfile('TMP','params'));
param0=param;
tic;
load(fullfile('TMP',sprintf('%d_params',nmod)),'param');
pscale=2^(iscale-1);
sbasis=param.basis;
if iscale>1
    sbasis='fem';
end
sizeim=ones(1,3);
switch sbasis
    case 'fem'

        mesh_file=fullfile('TMP',sprintf('%d_3d_mesh_%d',nmod,iscale-1));
        [dphidx,dphidy,dphidz]=CreateGradFiniteElementBasis3D(mesh_file,sizeim,pscale,[],'Gauss_points');
        wdetJ=GetWeigthDetJ(mesh_file,sizeim,1,'Gauss_points');

    case 'nurbs'
        mesh_file=fullfile('TMP',sprintf('%d_3d_mesh_%d',nmod,iscale-1));
        load(mesh_file,'degree');
            [dphidx,dphidy,dphidz,Xi,Yi,Zi,wdetJ]=CreateGradNURBSBasis3D(mesh_file,degree,'Gauss_points',1,'physical');

    otherwise
        error ('INVALID FUNCTIONAL BASIS');

end
        Nddl_tot=3*size(dphidx,2);
        phi0=sparse(size(dphidx,1),size(dphidx,2));
        epsxx=[dphidx,phi0,phi0];
        save(fullfile('TMP',sprintf('%d_3d_epsxx_%d',nmod,1*(iscale-1))),'epsxx','wdetJ','sizeim');
        epsyy=[phi0,dphidy,phi0];
        save(fullfile('TMP',sprintf('%d_3d_epsyy_%d',nmod,1*(iscale-1))),'epsyy','sizeim');
        epsxy=[dphidy,dphidx,phi0];
            Uxy=[dphidy,phi0,phi0];
            Uyx=[phi0,dphidx,phi0];
        save(fullfile('TMP',sprintf('%d_3d_epsxy_%d',nmod,1*(iscale-1))),'Uxy','Uyx','epsxy','sizeim');
        epsxz=[dphidz,phi0,dphidx];
        Uxz=[dphidz,phi0,phi0];
        Uzx=[phi0,phi0,dphidx];
        save(fullfile('TMP',sprintf('%d_3d_epsxz_%d',nmod,1*(iscale-1))),'Uxz','Uzx','epsxz','sizeim');
        epsyz=[phi0,dphidz,dphidy];
        Uyz=[phi0,dphidz,phi0];
        Uzy=[phi0,phi0,dphidy];
        save(fullfile('TMP',sprintf('%d_3d_epsyz_%d',nmod,1*(iscale-1))),'Uyz','Uzy','epsyz','sizeim');
        epszz=[phi0,phi0,dphidz];
        save(fullfile('TMP',sprintf('%d_3d_epszz_%d',nmod,1*(iscale-1))),'epszz','sizeim');

        if (iscale==1)&&(isfield(param,'enrichment'))

    mesh_file=fullfile('TMP',sprintf('%d_3d_mesh_%d',nmod,iscale-1));
    load(mesh_file,'ng','Nelems');
    wdetJo=GetWeigthDetJ(mesh_file,sizeim,1,'sub_cells',1);
    npt=size(wdetJo,1);
    epsxx=sparse(npt,size(epsxx,2));
    epsyy=sparse(npt,size(epsxx,2));
    epsxy=sparse(npt,size(epsxx,2));
    epsxz=sparse(npt,size(epsxx,2));
    epsyz=sparse(npt,size(epsxx,2));
    epszz=sparse(npt,size(epsxx,2));
    epscxx=sparse(npt,size(epsxx,2));
    epscyy=sparse(npt,size(epsxx,2));
    epscxy=sparse(npt,size(epsxx,2));
    epscxz=sparse(npt,size(epsxx,2));
    epscyz=sparse(npt,size(epsxx,2));
    epsczz=sparse(npt,size(epsxx,2));
        Uxy=sparse(npt,size(epsxx,2));
        Uyx=sparse(npt,size(epsxx,2));
        Uxz=sparse(npt,size(epsxx,2));
        Uzx=sparse(npt,size(epsxx,2));
        Uzy=sparse(npt,size(epsxx,2));
        Uyz=sparse(npt,size(epsxx,2));
        Ucxy=sparse(npt,size(epsxx,2));
        Ucyx=sparse(npt,size(epsxx,2));
        Ucxz=sparse(npt,size(epsxx,2));
        Uczx=sparse(npt,size(epsxx,2));
        Uczy=sparse(npt,size(epsxx,2));
        Ucyz=sparse(npt,size(epsxx,2));
    wdetJ=sparse(npt,1);
    if ~iscell(param.levelset_file)
        nbfis=1;
    else
        nbfis=numel(param.levelset_file);
    end
    roi=param0.roi;
    for ic=1:nbfis
        load(fullfile('TMP',sprintf('%d_enrichment_%d',nmod,ic)),'zone','face_nodes','face_elts','crackn','crackp');
        wdetJo=GetWeigthDetJ(mesh_file,sizeim,1,'sub_cells',face_elts);

        [dphidx,dphidy,dphidz]=CreateGradFiniteElementBasis3D(mesh_file,sizeim,pscale,face_nodes,'sub_cells',true,face_elts);
        phi0=sparse(size(dphidx,1),size(dphidx,2));
        epscxx=epscxx+[dphidx,phi0,phi0];
        epscyy=epscyy+[phi0,dphidy,phi0];
        epscxy=epscxy+[dphidy,dphidx,phi0];
        epscxz=epscxz+[dphidz,phi0,dphidx];
        epscyz=epscyz+[phi0,dphidz,dphidy];
        epsczz=epsczz+[phi0,phi0,dphidz];
        Ucxy=Ucxy+[dphidy,phi0,phi0];
        Ucyx=Ucyx+[phi0,dphidx,phi0];
        Ucxz=Ucxz+[dphidz,phi0,phi0];
        Uczx=Uczx+[phi0,phi0,dphidx];
        Uczy=Uczy+[phi0,phi0,dphidy];
        Ucyz=Ucyz+[phi0,dphidz,phi0];
        hn=double(crackn>=0);
        hn=diag(sparse(hn));
        heaviside=double(crackp>=0);
        heaviside=diag(sparse(heaviside));
        dphihdx=heaviside*dphidx(:,face_nodes)-dphidx(:,face_nodes)*hn;
        dphihdy=heaviside*dphidy(:,face_nodes)-dphidy(:,face_nodes)*hn;
        dphihdz=heaviside*dphidz(:,face_nodes)-dphidz(:,face_nodes)*hn;
        phih0=sparse(size(dphihdx,1),size(dphihdx,2));
        epsxx=[epsxx,dphihdx,phih0,phih0];
        epsyy=[epsyy,phih0,dphihdy,phih0];
        epsxy=[epsxy,dphihdy,dphihdx,phih0];
        epsxz=[epsxz,dphihdz,phih0,dphihdx];
        epsyz=[epsyz,phih0,dphihdz,dphihdy];
        epszz=[epszz,phih0,phih0,dphihdz];
        Uxy=[Uxy,dphihdy,phih0,phih0];
        Uyx=[Uyx,phih0,dphihdx,phih0];
        Uxz=[Uxz,dphihdz,phih0,phih0];
        Uzx=[Uzx,phih0,phih0,dphihdx];
        Uzy=[Uzy,phih0,phih0,dphihdy];
        Uyz=[Uyz,phih0,dphihdz,phih0];
        wdetJ=wdetJ+diag(wdetJo);
    end
        eps0=sparse(npt,size(epsxx,2)-size(epscxx,2));        
                   epscxx=[epscxx,eps0];
                    epscyy=[epscyy,eps0];
                   epscxy=[epscxy,eps0];
                   epscxz=[epscxz,eps0];
                   epscyz=[epscyz,eps0];
                   epsczz=[epsczz,eps0];
                   Ucxy=[Ucxy,eps0];
                   Ucxz=[Ucxz,eps0];
                   Uczy=[Uczy,eps0];
                   Ucyx=[Ucyx,eps0];
                   Uczx=[Uczx,eps0];
                   Ucyz=[Ucyz,eps0];
    wdetJ=diag(wdetJ);
    epsxx=epscxx+epsxx;
    epsyy=epscyy+epsyy;
    epsxy=epscxy+epsxy;
    epsxz=epscxz+epsxz;
    epsyz=epscyz+epsyz;
    epszz=epsczz+epszz;
    Uxy=Ucxy+Uxy;
    Uyx=Ucyx+Uyx;
    Uxz=Ucxz+Uxz;
    Uzx=Uczx+Uzx;
    Uzy=Uczy+Uzy;
    Uyz=Ucyz+Uyz;
    Nddl_tot=size(epsxx,2);
    save(fullfile('TMP',sprintf('%d_3d_xepsxx_%d',nmod,1*(iscale-1))),'epsxx','wdetJ','sizeim','Nddl_tot');
    save(fullfile('TMP',sprintf('%d_3d_xepsyy_%d',nmod,1*(iscale-1))),'epsyy','sizeim');
    save(fullfile('TMP',sprintf('%d_3d_xepsxy_%d',nmod,1*(iscale-1))),'epsxy','sizeim','Uxy','Uyx');
    save(fullfile('TMP',sprintf('%d_3d_xepsxz_%d',nmod,1*(iscale-1))),'epsxz','sizeim','Uxz','Uzx');
    save(fullfile('TMP',sprintf('%d_3d_xepsyz_%d',nmod,1*(iscale-1))),'epsyz','sizeim','Uzy','Uyz');
    save(fullfile('TMP',sprintf('%d_3d_xepszz_%d',nmod,1*(iscale-1))),'epszz','sizeim');
end
save(fullfile('TMP',sprintf('%d_3d_epsxx_%d',nmod,1*(iscale-1))),'Nddl_tot','-append');
disp(sprintf('Creating gradient basis function 3D for model %d...%6.2f s',nmod,toc));


end