source: issm/trunk/src/m/solvers/solver_stokescoupling_nonlinear.m@ 6581

function femmodel=solver_couplingstokes_nonlinear(femmodel,conserve_loads)
%SOLVER_COUPLINGSTOKES_NONLINEAR - core solver of coupling run
%
%   Usage:
%      [femmodel]=solver_couplingstokes_nonlinear(femmodel,conserve_loads)

        %Get parameters
        kffpartitioning=femmodel.parameters.Kff;

        %initialize solution vector
        converged=0; count=1;

        %First get ug=ug_horiz+ug_vert
        femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
        ug_horiz=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, femmodel.loads, femmodel.materials, femmodel.parameters);
        uf_horiz=Reducevectorgtof( ug_horiz, femmodel.nodesets,femmodel.parameters);

        while(~converged),

                %First compute the horizontal velocity
                femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);

                %Update the solution to make sure that vx and vxold are similar
                [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug_horiz);

                %save pointer to old velocity
                old_ug_horiz=ug_horiz;
                old_uf_horiz=uf_horiz;

                if kffpartitioning,
                        [K_gg_horiz,K_ff_horiz,K_fs_horiz,p_g_horiz,p_f_horiz,kmax_horiz]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
                        p_f_horiz = Reduceload( p_f_horiz, K_fs_horiz, femmodel.ys);
                else
                        [K_gg_horiz,K_ff_horiz,K_fs_horiz,p_g_horiz,p_f_horiz,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
                        [K_ff_horiz, K_fs_horiz] = Reducematrixfromgtof( K_gg_horiz, femmodel.nodesets); 
                        p_f_horiz = Reduceloadfromgtof( p_g_horiz, K_fs_horiz, femmodel.ys, femmodel.nodesets);
                end

                uf_horiz=Solver(K_ff_horiz,p_f_horiz,old_uf_horiz,femmodel.parameters);
                ug_horiz= Mergesolutionfromftog( uf_horiz, femmodel.ys, femmodel.nodesets,femmodel.parameters); 

                [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug_horiz);

                %Figure out if convergence have been reached
                converged=convergence(K_ff_horiz,p_f_horiz,uf_horiz,old_uf_horiz,femmodel.parameters);

                %Then compute vertical velocity
                femmodel=SetCurrentConfiguration(femmodel,DiagnosticVertAnalysisEnum);

                if kffpartitioning,
                        [K_gg_vert,K_ff_vert,K_fs_vert,p_g_vert,p_f_vert,kmax_vert]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
                        p_f_vert = Reduceload( p_f_vert, K_fs_vert, femmodel.ys);
                else
                        [K_gg_vert,K_ff_vert,K_fs_vert,p_g_vert,p_f_vert,kmax_vert]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
                        [K_ff_vert, K_fs_vert] = Reducematrixfromgtof( K_gg_vert, femmodel.nodesets); 
                        p_f_vert = Reduceloadfromgtof( p_g_vert, K_fs_vert, femmodel.ys, femmodel.nodesets);
                end

                uf_vert=Solver(K_ff_vert,p_f_vert,[],femmodel.parameters);
                ug_vert= Mergesolutionfromftog( uf_vert, femmodel.ys, femmodel.nodesets); 

                [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug_vert);

        end
end