Index: /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp (revision 4113) +++ /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp (revision 4114) @@ -10,5 +10,5 @@ #include "./solvers.h" -void solver_diagnostic_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, FemModel* fem,bool conserve_loads){ +void solver_diagnostic_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, FemModel* femmodel,bool conserve_loads){ @@ -42,17 +42,17 @@ /*Recover parameters: */ kflag=1; pflag=1; - fem->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); - fem->parameters->FindParam(&solver_string,SolverStringEnum); - fem->parameters->FindParam(&dim,DimEnum); - fem->parameters->FindParam(&verbose,VerboseEnum); - fem->parameters->FindParam(&min_mechanical_constraints,MinMechanicalConstraintsEnum); - fem->parameters->FindParam(&max_nonlinear_iterations,MaxNonlinearIterationsEnum); + femmodel->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); + femmodel->parameters->FindParam(&solver_string,SolverStringEnum); + femmodel->parameters->FindParam(&dim,DimEnum); + femmodel->parameters->FindParam(&verbose,VerboseEnum); + femmodel->parameters->FindParam(&min_mechanical_constraints,MinMechanicalConstraintsEnum); + femmodel->parameters->FindParam(&max_nonlinear_iterations,MaxNonlinearIterationsEnum); /*Were loads requested as output? : */ if(conserve_loads){ - loads=fem->loads->Copy(); //protect loads from being modified by the solution + loads=femmodel->loads->Copy(); //protect loads from being modified by the solution } else{ - loads=fem->loads; //modify loads in this solution + loads=femmodel->loads; //modify loads in this solution } @@ -61,6 +61,6 @@ /*Start non-linear iteration using input velocity: */ - GetSolutionFromInputsx(&ug, fem->elements, fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters); - Reducevectorgtofx(&uf, ug, fem->nodesets); + GetSolutionFromInputsx(&ug, femmodel->elements, femmodel->nodes, femmodel->vertices, loads, femmodel->materials, femmodel->parameters); + Reducevectorgtofx(&uf, ug, femmodel->nodesets); for(;;){ @@ -71,14 +71,11 @@ if (verbose) _printf_(" Generating matrices\n"); - //*Generate system matrices - SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag); + SystemMatricesx(&Kgg, &pg,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters,kflag,pflag); if (verbose) _printf_(" Generating penalty matrices\n"); - //*Generate penalty system matrices - PenaltySystemMatricesx(Kgg, pg,NULL,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag); - + PenaltySystemMatricesx(Kgg, pg,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters,kflag,pflag); + if (verbose) _printf_(" reducing matrix from g to f set\n"); - /*!Reduce matrix from g to f size:*/ - Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets); + Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->Gmn,femmodel->nodesets); /*Free ressources: */ @@ -86,6 +83,5 @@ if (verbose) _printf_(" reducing load from g to f set\n"); - /*!Reduce load from g to f size: */ - Reduceloadfromgtofx(&pf, pg, fem->Gmn, Kfs, fem->ys, fem->nodesets); + Reduceloadfromgtofx(&pf, pg, femmodel->Gmn, Kfs, femmodel->ys, femmodel->nodesets); //no need for pg and Kfs anymore @@ -93,27 +89,24 @@ MatFree(&Kfs); - /*Solve: */ if (verbose) _printf_(" solving\n"); Solverx(&uf, Kff, pf, old_uf, solver_string); - //Merge back to g set if (verbose) _printf_(" merging solution from f to g set\n"); - Mergesolutionfromftogx(&ug, uf,fem->Gmn,fem->ys,fem->nodesets); + Mergesolutionfromftogx(&ug, uf,femmodel->Gmn,femmodel->ys,femmodel->nodesets); //Update inputs using new solution: - InputUpdateFromSolutionx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,ug); + InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); - //Deal with penalty loads if (verbose) _printf_(" penalty constraints\n"); - PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters); + PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters); if(verbose)_printf_(" number of unstable constraints: %i\n",num_unstable_constraints); /*Figure out if convergence is reached.*/ - convergence(&converged,Kff,pf,uf,old_uf,fem->parameters); + convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters); MatFree(&Kff);VecFree(&pf); /*add converged to inputs: */ - InputUpdateFromConstantx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,converged,ConvergedEnum); + InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum); //rift convergence @@ -137,6 +130,6 @@ if (dim==3){ if(verbose)_printf_("%s\n","extruding velocity and pressure in 3d..."); - InputExtrudex( fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,VxEnum); - InputExtrudex( fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,VyEnum); + InputExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum); + InputExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum); } @@ -146,6 +139,6 @@ kflag=1; pflag=0; //stiffness generation only - SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag); - Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets); + SystemMatricesx(&Kgg, &pg,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters,kflag,pflag); + Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->Gmn,femmodel->nodesets); MatFree(&Kgg);VecFree(&pg); Index: /issm/trunk/src/m/solutions/jpl/NewFemModel.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/NewFemModel.m (revision 4113) +++ /issm/trunk/src/m/solutions/jpl/NewFemModel.m (revision 4114) @@ -31,17 +31,17 @@ displaystring(md.verbose,'%s',' generating single point constraints...'); - [femmodel.nodes,femmodel.yg(i)]=SpcNodes(femmodel.nodes,femmodel.constraints,analysis_type); + [femmodel.nodes,femmodel.m_yg(i)]=SpcNodes(femmodel.nodes,femmodel.constraints,analysis_type); displaystring(md.verbose,'%s',' generating rigid body constraints...'); - [femmodel.Rmg(i),femmodel.nodes]=MpcNodes(femmodel.nodes,femmodel.constraints,analysis_types); + [femmodel.m_Rmg(i),femmodel.nodes]=MpcNodes(femmodel.nodes,femmodel.constraints,analysis_types); displaystring(md.verbose,'%s',' generating node sets...'); - femmodel.nodesets(i)=BuildNodeSets(femmodel.nodes,analysis_type); + femmodel.m_nodesets(i)=BuildNodeSets(femmodel.nodes,analysis_type); displaystring(md.verbose,'%s',' reducing single point constraints vector...'); - femmodel.ys(i)=Reducevectorgtos(femmodel.yg(i).vector,femmodel.nodesets(i)); + femmodel.m_ys(i)=Reducevectorgtos(femmodel.m_yg(i).vector,femmodel.m_nodesets(i)); displaystring(md.verbose,'%s',' normalizing rigid body constraints matrix...'); - femmodel.Gmn(i)= NormalizeConstraints(femmodel.Rmg(i),femmodel.nodesets(i)); + femmodel.m_Gmn(i)= NormalizeConstraints(femmodel.m_Rmg(i),femmodel.m_nodesets(i)); displaystring(md.verbose,'%s',' configuring element and loads...'); Index: /issm/trunk/src/m/solutions/jpl/SetCurrentAnalysis.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/SetCurrentAnalysis.m (revision 4114) +++ /issm/trunk/src/m/solutions/jpl/SetCurrentAnalysis.m (revision 4114) @@ -0,0 +1,36 @@ +function femmodel=SetCurrentAnalysis(femmodel,analysis_enum) +%SETCURRENTANALYSIS - set current analysis used to configure elements and nodes in our solutions +% +% Usage: +% femmodel=SetCurrentAnalysis(femmodel,analysis_type) +% +% Ex: +% femmodel=SetCurrentAnalysis(femmodel,DiagnosticHorizAnalysisEnum) +% + + + %first, look for analysis: + int found=0; + for i=1:length(femmodel.analysis_type_list), + if analysis_type_list(i)==analysis_enum, + found=i; + break; + end + + if (found!=1), + error('SetCurrentAnalysis error message: could not find analysis_type in list of FemModel analyses'); + end + + + %activate matrices and vectors: + femmodel.Rmg=femmodel.m_Rmg(found); + femmodel.Gmn=femmodel.m_Gmn(found); + femmodel.nodesets=femmodel.m_nodesets(found); + femmodel.yg=femmodel.m_yg(found); + femmodel.ys=femmodel.m_ys(found); + + + %Now, plug analysis_counter and analysis_type inside the parameters: + %set counter and analyse_type + femmodel.parameters.AnalysisCounter=found; + femmodel.parameters.AnalysisType=analysis_enum; Index: /issm/trunk/src/m/solutions/jpl/SetCurrentAnalysisAlias.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/SetCurrentAnalysisAlias.m (revision 4114) +++ /issm/trunk/src/m/solutions/jpl/SetCurrentAnalysisAlias.m (revision 4114) @@ -0,0 +1,37 @@ +function femmodel=SetCurrentAnalysisAlias(femmodel,base_analysis_enum,real_analysis_type) +%SETCURRENTANALYSISALIAS- set current analysis used to configure elements and nodes in our solutions, +% and use another analysis_type instead. +% +% Usage: +% femmodel=SetCurrentAnalysisAlias(femmodel,base_analysis_type,real_analysis_type) +% +% Ex: +% femmodel=SetCurrentAnalysis(femmodel,DiagnosticHorizAnalysisEnum,AdjointAnalysisEnum) +% + + + %first, look for analysis: + int found=0; + for i=1:length(femmodel.analysis_type_list), + if analysis_type_list(i)==base_analysis_enum, + found=i; + break; + end + + if (found!=1), + error('SetCurrentAnalysisAlias error message: could not find analysis_type in list of FemModel analyses'); + end + + + %activate matrices and vectors: + femmodel.Rmg=femmodel.m_Rmg(found); + femmodel.Gmn=femmodel.m_Gmn(found); + femmodel.nodesets=femmodel.m_nodesets(found); + femmodel.yg=femmodel.m_yg(found); + femmodel.ys=femmodel.m_ys(found); + + + %Now, plug analysis_counter and analysis_type inside the parameters: + %set counter and analyse_type + femmodel.parameters.AnalysisCounter=found; + femmodel.parameters.AnalysisType=real_analysis_type; Index: sm/trunk/src/m/solutions/jpl/diagnostic_core_linear.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/diagnostic_core_linear.m (revision 4113) +++ (revision ) @@ -1,27 +1,0 @@ -function u_g=diagnostic_core_linear(m,analysis_type,sub_analysis_type) -%DIAGNOSTIC_CORE_LINEAR - linear solution sequence -% -% Usage: -% u_g=diagnostic_core_linear(m,analysis_type,sub_analysis_type) - - %stiffness and load generation only: - m.parameters.Kflag=1; m.parameters.Pflag=1; - - %system matrices - [K_gg, p_g]=SystemMatrices(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,analysis_type,sub_analysis_type); - [K_gg, p_g,kmax]=PenaltySystemMatrices(K_gg,p_g,m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,analysis_type,sub_analysis_type); - - %Reduce tangent matrix from g size to f size - [K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.Gmn, m.nodesets); - displaystring(m.parameters.Verbose>1,'%s%g',' condition number of stiffness matrix: ',condest(K_ff)); - - %Reduce load from g size to f size - [p_f] = Reduceloadfromgtof( p_g, m.Gmn, K_fs, m.ys, m.nodesets); - - %Solve - [u_f]=Solver(K_ff,p_f,[],m.parameters); - - %Merge back to g set - [u_g]= Mergesolutionfromftog( u_f, m.Gmn, m.ys, m.nodesets ); - -end %end function Index: sm/trunk/src/m/solutions/jpl/diagnostic_core_nonlinear.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/diagnostic_core_nonlinear.m (revision 4113) +++ (revision ) @@ -1,85 +1,0 @@ -function [u_g varargout]=diagnostic_core_nonlinear(m,input_loads,analysis_type,sub_analysis_type) -%DIAGNOSTIC_CORE_NONLINEAR - core solution of diagnostic non-linear -% -% Usage: -% [u_g varargout]=diagnostic_core_nonlinear(m,inputs_loads,analysis_type,sub_analysis_type) - -% stiffness and load generation only: - m.parameters.Kflag=1; m.parameters.Pflag=1; - -% Were loads requested as outputs ? - if isempty(input_loads), - loads=m.loads; - else - loads=input_loads; - end - -% initialize solution vector - converged=0; count=1; - - soln(count).u_g=GetSolutionFromInputs(m.elements, m.nodes, m.vertices, m.loads, m.materials, m.parameters, analysis_type, sub_analysis_type); - soln(count).u_f=Reducevectorgtof( soln(count).u_g, m.nodesets); - - displaystring(m.parameters.Verbose,'\n%s',[' starting direct shooting method']); - while(~converged), - - %system matrices - [K_gg_nopenalty , p_g_nopenalty]=SystemMatrices(m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,analysis_type,sub_analysis_type); - - %penalties - [K_gg , p_g, kmax]=PenaltySystemMatrices(K_gg_nopenalty,p_g_nopenalty,m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,analysis_type,sub_analysis_type); - - %Reduce tangent matrix from g size to f size - [K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.Gmn, m.nodesets); - - %Reduce load from g size to f size - [p_f] = Reduceloadfromgtof( p_g, m.Gmn, K_fs, m.ys, m.nodesets); - - %Increment index - count=count+1; - - %Solve - [soln(count).u_f]=Solver(K_ff,p_f,[],m.parameters); - - %Merge back to g set - [soln(count).u_g]= Mergesolutionfromftog( soln(count).u_f, m.Gmn, m.ys, m.nodesets ); - - %Update elements with new solution - m.elements=UpdateInputsFromSolution(m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,soln(count).u_g,analysis_type,sub_analysis_type); - - %penalty constraints - [loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints( m.elements,m.nodes,m.vertices,loads, m.materials,m.parameters,analysis_type,sub_analysis_type); - - displaystring(m.parameters.Verbose,'%s%i',' number of unstable constraints: ',num_unstable_constraints); - - %Figure out if convergence have been reached - converged=convergence(K_ff,p_f,soln(count).u_f,soln(count-1).u_f,m.parameters); - - %add convergence status into - [m.elements m.loads]=UpdateInputsFromConstant(m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,converged,ConvergedEnum); - - %rift convergence criterion - if ~constraints_converged, - converged=0; - end - - end - - %some cleanup - u_g=soln(count).u_g; - clear soln - - %more output might be needed, when running in control_core.m - nout=max(nargout,1)-1; - if nout==2, - inputs=add(inputs,'velocity',u_g,'doublevec',m.parameters.NumberOfDofsPerNode,m.parameters.NumberOfNodes); - m.parameters.Kflag=1; m.parameters.Pflag=0; - [K_gg, p_g]=SystemMatrices(m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,analysis_type,sub_analysis_type); - [K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.Gmn, m.nodesets); - varargout(1)={K_ff}; - varargout(2)={K_fs}; - end - if nout==1, - varargout(1)={loads}; - end -end Index: sm/trunk/src/m/solutions/jpl/thermal_core_nonlinear.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/thermal_core_nonlinear.m (revision 4113) +++ (revision ) @@ -1,74 +1,0 @@ -function [t_g ,loads, melting_offset]=thermal_core_nonlinear(m,analysis_type,sub_analysis_type) -%THERMAL_CORE_NONLINEAR - core of thermal solution sequence. -% model is return together with temperature -% -% Usage: -% [t_g ,loads, melting_offset]=thermal_core_nonlinear(m,analysis_type,sub_analysis_type); - - count=1; - converged=0; - -% we are going to return the loads, make them a variable of this routine - loads=m.loads; - - %stiffness and load generation only: - m.parameters.Kflag=1; m.parameters.Pflag=1; - - displaystring(m.parameters.Verbose,'\n%s',[' starting direct shooting method']); - - while(~converged), - - %Update inputs in datasets - if count==1 reset_penalties=0; else reset_penalties=1; end - [m.elements m.loads]=UpdateInputsFromConstant(m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,reset_penalties,ResetPenaltiesEnum); - - %system matrices - if ~m.parameters.Lowmem - if count==1 - displaystring(m.parameters.Verbose,'%s',[' system matrices']); - [K_gg_nopenalty, p_g_nopenalty]=SystemMatrices(m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,analysis_type,sub_analysis_type); - end - displaystring(m.parameters.Verbose,'%s',[' penalty system matrices']); - [K_gg , p_g, melting_offset]=PenaltySystemMatrices(K_gg_nopenalty,p_g_nopenalty,m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,analysis_type,sub_analysis_type); - else - m.parameters.Verbose - displaystring(m.parameters.Verbose,'%s',[' system matrices']); - [K_gg , p_g]=SystemMatrices(m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,analysis_type,sub_analysis_type); - displaystring(m.parameters.Verbose,'%s',[' penalty system matrices']); - [K_gg , p_g, melting_offset]=PenaltySystemMatrices(K_gg,p_g,m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,analysis_type,sub_analysis_type); - end - - %Reduce tangent matrix from g size to f size - [K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.Gmn, m.nodesets); - - %Reduce load from g size to f size - [p_f] = Reduceloadfromgtof( p_g, m.Gmn, K_fs, m.ys, m.nodesets); - - %Solve - displaystring(m.parameters.Verbose,'%s%g',' condition number of stiffness matrix: ',condest(K_ff)); - [t_f]=Solver(K_ff,p_f,[],m.parameters); - - %Merge back to g set - displaystring(m.parameters.Verbose,'%s',[' merging solution back to g set']); - [t_g]= Mergesolutionfromftog( t_f, m.Gmn, m.ys, m.nodesets ); - - %Update inputs in datasets - displaystring(m.parameters.Verbose,'%s',[' update inputs']); - [m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters] = UpdateInputsFromVector(m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,t_g,TemperatureEnum, VertexEnum); - - %penalty constraints - displaystring(m.parameters.Verbose,'%s',[' penalty constraints']); - [loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints(m.elements,m.nodes,m.vertices,loads, m.materials,m.parameters,analysis_type,sub_analysis_type); - - if ~converged, - displaystring(m.parameters.Verbose,'%s%i',' #unstable constraints ',num_unstable_constraints); - - if num_unstable_constraints<=m.parameters.MinThermalConstraints, - converged=1; - end - end - - count=count+1; - end - -end Index: /issm/trunk/src/m/solvers/solver_diagnostic_nonlinear.m =================================================================== --- /issm/trunk/src/m/solvers/solver_diagnostic_nonlinear.m (revision 4114) +++ /issm/trunk/src/m/solvers/solver_diagnostic_nonlinear.m (revision 4114) @@ -0,0 +1,103 @@ +function [femmodel ug varargout]=diagnostic_core_nonlinear(femmodel,conserve_loads) +%DIAGNOSTIC_CORE_NONLINEAR - core solution of diagnostic non-linear +% +% Usage: +% [femmodel ug varargout]=diagnostic_core_nonlinear(femmodel,conserve_loads) + + + %stiffness and load generation only: + femmodel.parameters.Kflag=1; femmodel.parameters.Pflag=1; + + %keep a copy of loads for now + loads=femmodel.loads; + + %initialize solution vector + converged=0; count=1; + + %Start non-linear iteration using input velocity: + ug=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, loads, femmodel.materials, femmodel.parameters); + uf=Reducevectorgtof( ug, femmodel.nodesets); + + while(~converged), + + %save pointer to old velocity + old_ug=ug; + old_uf=uf; + + displaystring(femmodel.parameters.Verbose,'%s',' Generating matrices:'); + [K_gg_nopenalty , p_g_nopenalty]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters); + + displaystring(femmodel.parameters.Verbose,'%s',' Generating penalty matrices:'); + [K_gg , p_g, kmax]=PenaltySystemMatrices(K_gg_nopenalty,p_g_nopenalty,femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters); + + displaystring(femmodel.parameters.Verbose,'%s',' reducing matrix from g to f set'); + [K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.Gmn, femmodel.nodesets); + + displaystring(femmodel.parameters.Verbose,'%s',' reduce load from g size to f size'); + p_f = Reduceloadfromgtof( p_g, femmodel.Gmn, K_fs, femmodel.ys, femmodel.nodesets); + + displaystring(femmodel.parameters.Verbose,'%s',' solving'); + uf=Solver(K_ff,p_f,old_uf,femmodel.parameters); + + displaystring(femmodel.parameters.Verbose,'%s',' merge back to g set'); + ug= Mergesolutionfromftog( uf, femmodel.Gmn, femmodel.ys, femmodel.nodesets ); + + %Update inputs using new solution + femmodel.elements=UpdateInputsFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,ug); + + displaystring(femmodel.parameters.Verbose,'%s',' penalty constraints'); + [loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints( femmodel.elements,femmodel.nodes,femmodel.vertices,loads, femmodel.materials,femmodel.parameters); + + displaystring(femmodel.parameters.Verbose,'%s%i',' number of unstable constraints: ',num_unstable_constraints); + + %Figure out if convergence have been reached + converged=convergence(K_ff,p_f,uf,old_uf,soln(count-1).uf,femmodel.parameters); + + %add convergence status into status + [femmodel.elements loads]=UpdateInputsFromConstant(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,converged,ConvergedEnum); + + %rift convergence + if ~constraints_converged, + if converged, + if num_unstable_constraints <= min_mechanical_constraints, + converged=1; + end + else + converged=0; + end + end + end + + %increase count + count=count+1; + if(converged==1)break; + if(count>max_nonlinear_iterations), + displaystring(femmodel.parameters.Verbose,'%s%i%s',' maximum number of iterations ',max_nonlinear_iterations,' exceeded'); + end + end + + end + + %extrude if we are in 3D: + if dim==3, + displaystring(femmodel.parameters.Verbose,'%s',' extruding velocity and pressure in 3d...'); + InputExtrude( femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VxEnum); + InputExtrude( femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VyEnum); + end + + %more output might be needed, when running in control_core.m + nout=max(nargout,1)-2; + if nout==2, + femmodel.parameters.Kflag=1; femmodel.parameters.Pflag=0; + [K_gg, p_g]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters); + [K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.Gmn, femmodel.nodesets); + varargout(1)={K_ff}; + varargout(2)={K_fs}; + end + + + %deal with loads: + if conserve_loads==false, + femmodel.loads=loads; + end +end Index: /issm/trunk/src/m/solvers/solver_linear.m =================================================================== --- /issm/trunk/src/m/solvers/solver_linear.m (revision 4114) +++ /issm/trunk/src/m/solvers/solver_linear.m (revision 4114) @@ -0,0 +1,27 @@ +function u_g=diagnostic_core_linear(m,analysis_type,sub_analysis_type) +%DIAGNOSTIC_CORE_LINEAR - linear solution sequence +% +% Usage: +% u_g=diagnostic_core_linear(m,analysis_type,sub_analysis_type) + + %stiffness and load generation only: + m.parameters.Kflag=1; m.parameters.Pflag=1; + + %system matrices + [K_gg, p_g]=SystemMatrices(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,analysis_type,sub_analysis_type); + [K_gg, p_g,kmax]=PenaltySystemMatrices(K_gg,p_g,m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,analysis_type,sub_analysis_type); + + %Reduce tangent matrix from g size to f size + [K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.Gmn, m.nodesets); + displaystring(m.parameters.Verbose>1,'%s%g',' condition number of stiffness matrix: ',condest(K_ff)); + + %Reduce load from g size to f size + [p_f] = Reduceloadfromgtof( p_g, m.Gmn, K_fs, m.ys, m.nodesets); + + %Solve + [u_f]=Solver(K_ff,p_f,[],m.parameters); + + %Merge back to g set + [u_g]= Mergesolutionfromftog( u_f, m.Gmn, m.ys, m.nodesets ); + +end %end function Index: /issm/trunk/src/m/solvers/solver_thermal_nonlinear.m =================================================================== --- /issm/trunk/src/m/solvers/solver_thermal_nonlinear.m (revision 4114) +++ /issm/trunk/src/m/solvers/solver_thermal_nonlinear.m (revision 4114) @@ -0,0 +1,74 @@ +function [t_g ,loads, melting_offset]=thermal_core_nonlinear(m,analysis_type,sub_analysis_type) +%THERMAL_CORE_NONLINEAR - core of thermal solution sequence. +% model is return together with temperature +% +% Usage: +% [t_g ,loads, melting_offset]=thermal_core_nonlinear(m,analysis_type,sub_analysis_type); + + count=1; + converged=0; + +% we are going to return the loads, make them a variable of this routine + loads=m.loads; + + %stiffness and load generation only: + m.parameters.Kflag=1; m.parameters.Pflag=1; + + displaystring(m.parameters.Verbose,'\n%s',[' starting direct shooting method']); + + while(~converged), + + %Update inputs in datasets + if count==1 reset_penalties=0; else reset_penalties=1; end + [m.elements m.loads]=UpdateInputsFromConstant(m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,reset_penalties,ResetPenaltiesEnum); + + %system matrices + if ~m.parameters.Lowmem + if count==1 + displaystring(m.parameters.Verbose,'%s',[' system matrices']); + [K_gg_nopenalty, p_g_nopenalty]=SystemMatrices(m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,analysis_type,sub_analysis_type); + end + displaystring(m.parameters.Verbose,'%s',[' penalty system matrices']); + [K_gg , p_g, melting_offset]=PenaltySystemMatrices(K_gg_nopenalty,p_g_nopenalty,m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,analysis_type,sub_analysis_type); + else + m.parameters.Verbose + displaystring(m.parameters.Verbose,'%s',[' system matrices']); + [K_gg , p_g]=SystemMatrices(m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,analysis_type,sub_analysis_type); + displaystring(m.parameters.Verbose,'%s',[' penalty system matrices']); + [K_gg , p_g, melting_offset]=PenaltySystemMatrices(K_gg,p_g,m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,analysis_type,sub_analysis_type); + end + + %Reduce tangent matrix from g size to f size + [K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.Gmn, m.nodesets); + + %Reduce load from g size to f size + [p_f] = Reduceloadfromgtof( p_g, m.Gmn, K_fs, m.ys, m.nodesets); + + %Solve + displaystring(m.parameters.Verbose,'%s%g',' condition number of stiffness matrix: ',condest(K_ff)); + [t_f]=Solver(K_ff,p_f,[],m.parameters); + + %Merge back to g set + displaystring(m.parameters.Verbose,'%s',[' merging solution back to g set']); + [t_g]= Mergesolutionfromftog( t_f, m.Gmn, m.ys, m.nodesets ); + + %Update inputs in datasets + displaystring(m.parameters.Verbose,'%s',[' update inputs']); + [m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters] = UpdateInputsFromVector(m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,t_g,TemperatureEnum, VertexEnum); + + %penalty constraints + displaystring(m.parameters.Verbose,'%s',[' penalty constraints']); + [loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints(m.elements,m.nodes,m.vertices,loads, m.materials,m.parameters,analysis_type,sub_analysis_type); + + if ~converged, + displaystring(m.parameters.Verbose,'%s%i',' #unstable constraints ',num_unstable_constraints); + + if num_unstable_constraints<=m.parameters.MinThermalConstraints, + converged=1; + end + end + + count=count+1; + end + +end