function run_fem_job(param,model,Up,nmod)
if nargin<3,Up=[];end
if nargin<3,nmod=0;end
clear functions
pshift=0;
if isfield(param,'deformed_image')
if iscell(param.deformed_image)
    if size(param.deformed_image,1)>2
        pshift=size(param.deformed_image,1);
        defs=param.deformed_image;
        param.deformed_image=defs(1,:);
        if size(defs,2)==1,param.deformed_image=defs{1,1};end
    end
end
end
param.onflight=1;
LoadParameters(param);
LoadParameters(model,nmod);
ReferenceImage(nmod);
LoadMask(nmod);
LoadMeshes(nmod);
LoadMat(nmod);

nscale=model.nscale;
filres=param.result_file;
%%
load(fullfile('TMP','sample0_0.mat'),'sizeim')
if prod(sizeim)>4e6
    param.initialization='bilin+';
    %    param.initialization='rbt_elem';
    LoadParameters(param);
end
load(fullfile('TMP','sample0_0.mat'),'sizeim')
% if isfield(param,'detect')
%     if param.detect
%      param.initialization='rbt';
%     LoadParameters(param);
%     end
% end
load(fullfile('TMP',sprintf('%d_mesh_%d',nmod,1-1)),'selected','Nnodes')
if isfield(model,'mesh_type')
    if model.mesh_type==1
        rflag=true;
        save(fullfile('TMP',sprintf('%d_mesh_%d',nmod,1-1)),'rflag','-append')
    end
end
extract=0;
if strcmp(param.regularization_type,'equilibrium_gap')
    for iz=1:size(model.zone,2)
        zone=model.zone(:,iz);
        switch zone{4}
            %            case 6
            %                selected(zone{6})=0;
            case 5
                if (zone{8}>0)
                    extract=1;
                end
        end
    end
end
tips=zeros(size(model.zone,2),2);
cracks=[];
inds=[];
if extract
    load(fullfile('TMP',sprintf('%d_mesh_%d',nmod,1-1)),'xo','yo','conn','elt')
    if param.detect
        modes=1:2;
        harms=-3:7;
        
    else
        modes=1:2;
        harms=0:7;
    end
    roi=param.roi;
    freenodes=ones(length(xo),1);
    xo=xo+roi(1)-1;
    yo=yo+roi(3)-1;
    L=zeros(prod(Nnodes),1);
    nu=model.material_parameters.nu;
    kappa=(3-4*nu);
    %    kappa=(3-nu)/(1+nu)
    
    zo=xo+1i*yo;
    nw=0;
    ntip=0;
    for iz=1:size(model.zone,2)
        zone=model.zone(:,iz);
        if zone{4}==5
            if (zone{8}>0)
                tips(iz,1)=ntip+1;
                indc=zone{10};
                cnodes=indc{1};
                ztip=zo(cnodes(1,1));
                tt=-diff(zo(cnodes(1:2,1)));
                tt=tt/abs(tt);
                rtip=mean(abs(zo(indc{2})-ztip));
                
                %                rtip=Inf
                
                in=abs(zo-ztip)<=rtip;
                in(indc{2})=1;
                freenodes(in)=0;
                
                
                selected(in)=0;
                selected(indc{2})=1;
                
                
                zzone=(zo(in)-ztip)*tt';
                idin=1:prod(Nnodes);
                idin=idin(in);
                for ip=1:length(zzone)
                    if any(cnodes(:,1)==idin(ip))
                        zzone(ip)=abs(zzone(ip))*exp(1i*pi);
                    end
                    if any(cnodes(:,2)==idin(ip))
                        zzone(ip)=abs(zzone(ip))*exp(-1i*pi);
                    end
                end
                phi=Williams(zzone,modes,harms,kappa);
                phi=phi*tt;
                
                L(in,nw+(1:size(phi,2)))=phi;
                inds(nw+(1:length(harms)))=harms;
                ntip=ntip+1;
                nw=nw+size(phi,2);
                if  (zone{9}>0)
                    tips(iz,2)=ntip+1;
                    ztip=zo(cnodes(end,1));
                    tt=diff(zo(cnodes(end+(-1:0),1)));
                    tt=tt/abs(tt);
                    rtip=mean(abs(zo(indc{3})-ztip));
                    in=abs(zo-ztip)<=rtip;
                    in(indc{3})=1;
                    freenodes(in)=0;
                    
                    selected(in)=0;
                    selected(indc{3})=1;
                    
                    zzone=(zo(in)-ztip)*tt';
                    idin=1:prod(Nnodes);
                    idin=idin(in);
                    for ip=1:length(zzone)
                        if any(cnodes(:,1)==idin(ip))
                            zzone(ip)=abs(zzone(ip))*exp(-1i*pi);
                        end
                        if any(cnodes(:,2)==idin(ip))
                            zzone(ip)=abs(zzone(ip))*exp(1i*pi);
                        end
                    end
                    phi=Williams(zzone,modes,harms,kappa);
                    phi=phi*tt;
                    L(in,nw+(1:size(phi,2)))=phi;
                    inds(nw+(1:length(harms)))=harms;
                    nw=nw+size(phi,2);
                    ntip=ntip+1;
                end
                
            end
        end
    end
    nfree=sum(freenodes);
    indK=2*nfree+(1:size(L,2));
    Ln=sparse(find(freenodes),1:nfree,1,prod(Nnodes),nfree);
    L=[real(L);imag(L)];
    L=[blkdiag(Ln,Ln),L];
    %    selected(~freenodes)=0;
    
    if param.detect
        tnodes=[];
        for iz=1:size(model.zone,2)
            zone=model.zone(:,iz);
            if zone{4}==5
                if (zone{8}>0)
                    indc=zone{10};
                    cnodes=indc{1};
                    %                     htip=zone{7};
                    %                     rtip=zone{8};
                    
                    ztip=zo(cnodes(1,1));
                    rtip=min(abs(zo(indc{4})-ztip));
                    tnodes=abs(zo-ztip)<rtip;
                    
                    
                    %                     ntip=ceil(min(20,0.2*rtip)/htip);
                    %                     ttnodes=cnodes(1);
                    %                     for il=1:ntip
                    %                         ttnodes=AddOneNodeLayer(conn,ttnodes);
                    %                     end
                    %                         tnodes=[tnodes;ttnodes];
                    if  (zone{9}>0)
                        ztip=zo(cnodes(end,1));
                        rtip=min(abs(zo(indc{5})-ztip));
                        tnodes=tnodes|(abs(zo-ztip)<rtip);
                        %                                     rtip=mean(abs(zo(indc{3})-ztip));
                        %                 ttnodes=abs(zo-ztip)<=0.4*rtip;
                        
                        
                        %                         ttnodes=cnodes(end,1);
                        %                         for il=1:ntip
                        %                             ttnodes=AddOneNodeLayer(conn,ttnodes);
                        %                         end
                        %                         tnodes=[tnodes;ttnodes];
                    end
                end
            end
        end
        
        ielt=GetEltsFromNodes(conn,elt,tnodes,0);
        elt(ielt)=[];
        conn=conn(:,1:3);
        conn(ielt,:)=[];
        keep=unique(conn(:));
        xo=xo(keep)-roi(1)+1;
        yo=yo(keep)-roi(3)+1;
        selected=selected(keep);
        L=L([keep(:);keep+prod(Nnodes)],:);
        
        newids=zeros(prod(Nnodes),1);
        newids(keep)=1:length(keep);
        conn=newids(conn(:,1:3));
        
        elt=3*ones(size(conn,1),1);
        Nnodes=[length(xo),1,1];
        Nelems=[length(elt),1,1];
        conn=[conn,zeros(length(elt),1)];
        save(fullfile('TMP',sprintf('%d_mesh_%d',nmod,1-1)),'xo','yo','conn','elt','Nnodes','Nelems','-append')
    end
end
% param.regularization_type='none';
% param0.regularization_parameter=64;
% LoadParameters(param)
save(fullfile('TMP',sprintf('%d_mesh_%d',nmod,1-1)),'selected','-append')
%%
U1=[];Ks=[];
if ~isempty(Up),nscale=1;end
if param.thermo==1
    if pshift, error('NOT CODED YET');end
    for iscale=nscale:-1:1
        if isempty(Up)
            Uini=InitializeSolution(U1,iscale,nmod);
        else
            Uini=Up;
        end
        
        if (iscale==1)&&extract
            Pini=(L'*L)\(L'*Uini);
            %            Pini(inds<0,:)=0;
            Uini=L*Pini;
            [U1,dKs]=SolveIR2D(Uini,iscale,nmod,L);
            Ks=dKs(indK,:);
        else
            [U1]=SolveIR2D(Uini,iscale,nmod);
        end
    end
else
    if prod(sizeim)<4e6&&~param.normalize_grey_level&&(param.psample==1)
        
        for iscale=nscale:-1:1
            disp(sprintf('Pre-processing scale %d...',iscale));
            CreateBasisFunction(iscale,nmod);
            ComputeGradFPhi(iscale,nmod);
            CreateGradBasisFunction(iscale,nmod);
            AssembleCorrelationOperator(iscale,nmod);
            
            %%
            
            
            if isempty(Up)
                Uini=InitializeSolution(U1,iscale,nmod);
            else
                Uini=Up;
            end
            if (iscale==1)&&extract
                if pshift, error('NOT CODED YET');end
                Pini=(L'*L)\(L'*Uini);
                %            Pini(indK(inds<0),:)=0;
                %Pini(indK)'
                %keyboard
                %pause
                Uini=L*Pini;
                [U1,dKs]=RBMSolver(Uini,iscale,nmod,L);
                Ks=Pini(indK,:)+dKs(indK,:);
            else
                [U1]=Solve(Uini,iscale,nmod);
                
                if pshift&&iscale==1
                    Uini=U1;
                    for icam=2:pshift
                        if size(defs,2)==1
                            param.deformed_image=defs{icam,1};
                        else
                            param.deformed_image=defs(icam,:);
                        end
                        param.restart=1;
                        LoadParameters(param);
                        [Ui]=Solve(Uini,iscale,nmod);
                        Ux=mean(Uini((1:prod(Nnodes)))-Ui((1:prod(Nnodes))));
                        Uy=mean(Uini(prod(Nnodes)+(1:prod(Nnodes)))-Ui(prod(Nnodes)+(1:prod(Nnodes))));
                        Ui((1:prod(Nnodes)))=Ui((1:prod(Nnodes)))+Ux;
                        Ui(prod(Nnodes)+(1:prod(Nnodes)))=Ui(prod(Nnodes)+(1:prod(Nnodes)))+Uy;
                        U1=U1+Ui;
                    end
                    U1=U1/pshift;
                    param.iter_max=2;
                    if size(defs,2)==1
                        param.deformed_image=defs{1,1};
                    else
                        param.deformed_image=defs(1,:);
                    end
                    LoadParameters(param)
                    [~]=Solve(U1,1,nmod);
                    
                end
            end
            
        end
        
    else
        for iscale=nscale:-1:1
            if isempty(Up)
                Uini=InitializeSolution(U1,iscale,nmod);
            else
                Uini=Up;
            end
            
            if (iscale==1)&&extract
                if pshift, error('NOT CODED YET');end
                Pini=(L'*L)\(L'*Uini);
                %            Pini(inds<0,:)=0;
                Uini=L*Pini;
                [U1,dKs]=Solve2D(Uini,iscale,nmod,L);
                Ks=dKs(indK,:);
            else
                [U1]=Solve2D(Uini,iscale,nmod);
                if pshift&&iscale==1
                    Uini=U1;
                    for icam=2:pshift
                        if size(defs,2)==1
                            param.deformed_image=defs{icam,1};
                        else
                            param.deformed_image=defs(icam,:);
                        end
                        param.restart=1;
                        LoadParameters(param);
                        [Ui]=Solve2D(Uini,iscale,nmod);
                        Ux=mean(Uini((1:prod(Nnodes)))-Ui((1:prod(Nnodes))));
                        Uy=mean(Uini(prod(Nnodes)+(1:prod(Nnodes)))-Ui(prod(Nnodes)+(1:prod(Nnodes))));
                        Ui((1:prod(Nnodes)))=Ui((1:prod(Nnodes)))+Ux;
                        Ui(prod(Nnodes)+(1:prod(Nnodes)))=Ui(prod(Nnodes)+(1:prod(Nnodes)))+Uy;
                        U1=U1+Ui;
                    end
                    U1=U1/pshift;
                    param.iter_max=2;
                    if size(defs,2)==1
                        param.deformed_image=defs{1,1};
                    else
                        param.deformed_image=defs(1,:);
                    end
                    LoadParameters(param)
                    [~]=Solve2D(U1,1,nmod);
                    
                end
                
            end
        end
    end
end
%%
filreso=strrep(filres,'.res','');
U=U1;
copyfile(fullfile('TMP',sprintf('%d_error_%d.mat',nmod,1-1)), [filreso,'-error.res']);
load(fullfile('TMP',sprintf('%d_params',nmod)),'param');
model=param;
load(fullfile('TMP','params'),'param');
load(fullfile('TMP',sprintf('%d_mesh_%d',nmod,1-1)),'Nnodes','Nelems','xo','yo','conn','elt','ng','rflag','rint');
try
    tips(~(cell2mat(model.zone(4,:))==5),:)=[];
    model.zone(:,~(cell2mat(model.zone(4,:))==5))=[];
    cracks=find(cell2mat(model.zone(4,:))==5);
catch, end
%ExportImageToVTK(filres)
delete([fullfile('VTK',[filreso,'-error']),'-0*.vtk']);
if prod(sizeim)<4e6&~param.normalize_grey_level&&(param.thermo==0)&&(param.psample==1)
    if isfield(param,'deformed_image')
        if isfield(param,'image_number')
            images=param.image_number;
        else
            if iscell(param.deformed_image)
                images=1:size(param.deformed_image,2);
            else
                images=1;
            end
        end
    else
        images=(2:param.video_sampling:param.number_of_frames)-1;
    end
    for iim=1:size(U,2)
        movefile(fullfile('VTK',sprintf('camr-1-camd-1-scale-1-%d-error.vtk',iim)),fullfile('VTK',sprintf('%s-error-%04d.vtk',filreso,images(iim))));
    end
    VTKExportScalarMap([filreso,sprintf('-%d',0)],'error',param.roi(1)-1+(1:sizeim(1)),param.roi(3)-1+(1:sizeim(2)),zeros(sizeim),1);
    movefile(fullfile('VTK',[filreso,sprintf('-%d',0),'-error.vtk']),fullfile('VTK',sprintf('%s-error-%04d.vtk',filreso,0)));
end
delete(fullfile('VTK',['camr*','-error.vtk']));
save(filres,'U','Ks','tips','cracks','Nnodes','Nelems','xo','yo','param','model','nmod','conn','elt','rint','ng','rflag','-v7.3');
if isfield(param,'calibration_data')
    [U,Xo,Yo,Zo]=Extract3DFields(U1,nmod);
    CreateGradBasisFunction25D(iscale,nmod);
    
    save(filres,'U','U1','Xo','Yo','Zo','-append');
    writeVTKmesh25D(filres);
    postproVTK25D(filres,1);
else
    postproVTK([filres,''],0);
end
end