Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 5996) +++ /issm/trunk/src/c/Makefile.am (revision 5997) @@ -1146,6 +1146,4 @@ ./solvers/solver_adjoint_linear.cpp\ ./solvers/solver_diagnostic_nonlinear.cpp\ - ./solvers/solver_diagnostic_nonlinear_kgg.cpp\ - ./solvers/solver_diagnostic_nonlinear_kff.cpp\ ./solvers/solver_thermal_nonlinear.cpp\ ./modules/Bamgx/Bamgx.cpp\ Index: /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp (revision 5996) +++ /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp (revision 5997) @@ -12,12 +12,91 @@ void solver_diagnostic_nonlinear(FemModel* femmodel,bool conserve_loads){ + /*intermediary: */ + Mat Kgg = NULL, Kff = NULL, Kfs = NULL; + Vec ug = NULL, uf = NULL, old_ug= NULL, old_uf = NULL; + Vec pg = NULL, pf = NULL; + Loads* loads=NULL; + int converged; + int constraints_converged; + int num_unstable_constraints; + int count; + /*parameters:*/ + int verbose=0; bool kff=false; + int min_mechanical_constraints; + int max_nonlinear_iterations; /*Recover parameters: */ + femmodel->parameters->FindParam(&verbose,VerboseEnum); femmodel->parameters->FindParam(&kff,KffEnum); + femmodel->parameters->FindParam(&min_mechanical_constraints,MinMechanicalConstraintsEnum); + femmodel->parameters->FindParam(&max_nonlinear_iterations,MaxNonlinearIterationsEnum); + + /*Were loads requested as output? : */ + if(conserve_loads) loads=(Loads*)femmodel->loads->Copy(); //protect loads from being modified by the solution + else loads=(Loads*)femmodel->loads; //modify loads in this solution - if(!kff) solver_diagnostic_nonlinear_kgg(femmodel,conserve_loads); - else solver_diagnostic_nonlinear_kff(femmodel,conserve_loads); + count=1; + converged=0; + /*Start non-linear iteration using input velocity: */ + GetSolutionFromInputsx(&ug, femmodel->elements, femmodel->nodes, femmodel->vertices, loads, femmodel->materials, femmodel->parameters); + Reducevectorgtofx(&uf, ug, femmodel->nodesets,femmodel->parameters); + + //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate) + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); + + for(;;){ + + //save pointer to old velocity + VecFree(&old_ug);old_ug=ug; + VecFree(&old_uf);old_uf=uf; + + if(kff){ + SystemMatricesx(NULL,&Kff, &Kfs, NULL,&pf, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters); + Reduceloadx(pf, Kfs, femmodel->ys, femmodel->parameters); MatFree(&Kfs); + } + else{ + SystemMatricesx(&Kgg, NULL, NULL, &pg,NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters); + Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets,femmodel->parameters); MatFree(&Kgg); + Reduceloadfromgtofx(&pf, pg, Kfs, femmodel->ys, femmodel->nodesets,femmodel->parameters); VecFree(&pg); MatFree(&Kfs); + } + + Solverx(&uf, Kff, pf, old_uf, femmodel->parameters); + + Mergesolutionfromftogx(&ug, uf,femmodel->ys,femmodel->nodesets,femmodel->parameters); + + InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); + + 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); + + convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters); MatFree(&Kff);VecFree(&pf); + + InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum); + + //rift convergence + if (!constraints_converged) { + if (converged){ + if (num_unstable_constraints <= min_mechanical_constraints) converged=1; + else converged=0; + } + } + + /*Increase count: */ + count++; + if(converged==1)break; + if(count>=max_nonlinear_iterations){ + _printf_(" maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations); + break; + } + } + + /*clean-up*/ + if(conserve_loads) delete loads; + VecFree(&uf); + VecFree(&ug); + VecFree(&old_uf); + VecFree(&old_ug); } Index: sm/trunk/src/c/solvers/solver_diagnostic_nonlinear_kff.cpp =================================================================== --- /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear_kff.cpp (revision 5996) +++ (revision ) @@ -1,99 +1,0 @@ -/*!\file: solver_diagnostic_nonlinear_kff.cpp - * \brief: core of the diagnostic solution for non linear materials, using f and s-set - */ - -#include "../toolkits/toolkits.h" -#include "../objects/objects.h" -#include "../EnumDefinitions/EnumDefinitions.h" -#include "../modules/modules.h" -#include "../solutions/solutions.h" -#include "./solvers.h" - -void solver_diagnostic_nonlinear_kff(FemModel* femmodel,bool conserve_loads){ - - /*intermediary: */ - Mat Kff = NULL, Kfs = NULL; - Vec ug = NULL, uf = NULL, old_ug= NULL, old_uf = NULL; - Vec pf = NULL; - Loads* loads=NULL; - int converged; - int constraints_converged; - int num_unstable_constraints; - int count; - int min_mechanical_constraints; - int max_nonlinear_iterations; - - /*parameters:*/ - int verbose=0; - - /*Recover parameters: */ - 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=(Loads*)femmodel->loads->Copy(); //protect loads from being modified by the solution - } - else{ - loads=(Loads*)femmodel->loads; //modify loads in this solution - } - - count=1; - converged=0; - - /*Start non-linear iteration using input velocity: */ - GetSolutionFromInputsx(&ug, femmodel->elements, femmodel->nodes, femmodel->vertices, loads, femmodel->materials, femmodel->parameters); - Reducevectorgtofx(&uf, ug, femmodel->nodesets,femmodel->parameters); - - //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate) - InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); - - for(;;){ - - //save pointer to old velocity - VecFree(&old_ug);old_ug=ug; - VecFree(&old_uf);old_uf=uf; - - SystemMatricesx(NULL,&Kff, &Kfs, NULL,&pf, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters); - - Reduceloadx(pf, Kfs, femmodel->ys, femmodel->parameters); MatFree(&Kfs); - - Solverx(&uf, Kff, pf, old_uf, femmodel->parameters); - - Mergesolutionfromftogx(&ug, uf,femmodel->ys,femmodel->nodesets,femmodel->parameters); - - InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); - - 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); - - convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters); MatFree(&Kff);VecFree(&pf); - - InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum); - - //rift convergence - if (!constraints_converged) { - if (converged){ - if (num_unstable_constraints <= min_mechanical_constraints) converged=1; - else converged=0; - } - } - - /*Increase count: */ - count++; - if(converged==1)break; - if(count>=max_nonlinear_iterations){ - _printf_(" maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations); - break; - } - } - - /*clean-up*/ - if(conserve_loads) delete loads; - VecFree(&uf); - VecFree(&ug); - VecFree(&old_uf); - VecFree(&old_ug); - -} Index: sm/trunk/src/c/solvers/solver_diagnostic_nonlinear_kgg.cpp =================================================================== --- /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear_kgg.cpp (revision 5996) +++ (revision ) @@ -1,101 +1,0 @@ -/*!\file: solver_diagnostic_nonlinear.cpp - * \brief: core of the diagnostic solution for non linear materials - */ - -#include "../toolkits/toolkits.h" -#include "../objects/objects.h" -#include "../EnumDefinitions/EnumDefinitions.h" -#include "../modules/modules.h" -#include "../solutions/solutions.h" -#include "./solvers.h" - -void solver_diagnostic_nonlinear_kgg(FemModel* femmodel,bool conserve_loads){ - - /*intermediary: */ - Mat Kgg = NULL, Kff = NULL, Kfs = NULL; - Vec ug = NULL, uf = NULL, old_ug= NULL, old_uf = NULL; - Vec pg = NULL, pf = NULL; - Loads* loads=NULL; - int converged; - int constraints_converged; - int num_unstable_constraints; - int count; - int min_mechanical_constraints; - int max_nonlinear_iterations; - - /*parameters:*/ - int verbose=0; - - /*Recover parameters: */ - 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=(Loads*)femmodel->loads->Copy(); //protect loads from being modified by the solution - } - else{ - loads=(Loads*)femmodel->loads; //modify loads in this solution - } - - count=1; - converged=0; - - /*Start non-linear iteration using input velocity: */ - GetSolutionFromInputsx(&ug, femmodel->elements, femmodel->nodes, femmodel->vertices, loads, femmodel->materials, femmodel->parameters); - Reducevectorgtofx(&uf, ug, femmodel->nodesets,femmodel->parameters); - - //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate) - InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); - - for(;;){ - - //save pointer to old velocity - VecFree(&old_ug);old_ug=ug; - VecFree(&old_uf);old_uf=uf; - - SystemMatricesx(&Kgg, NULL, NULL, &pg,NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters); - - Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets,femmodel->parameters); MatFree(&Kgg); - - Reduceloadfromgtofx(&pf, pg, Kfs, femmodel->ys, femmodel->nodesets,femmodel->parameters); VecFree(&pg); MatFree(&Kfs); - - Solverx(&uf, Kff, pf, old_uf, femmodel->parameters); - - Mergesolutionfromftogx(&ug, uf,femmodel->ys,femmodel->nodesets,femmodel->parameters); - - InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); - - 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); - - convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters); MatFree(&Kff);VecFree(&pf); - - InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum); - - //rift convergence - if (!constraints_converged) { - if (converged){ - if (num_unstable_constraints <= min_mechanical_constraints) converged=1; - else converged=0; - } - } - - /*Increase count: */ - count++; - if(converged==1)break; - if(count>=max_nonlinear_iterations){ - _printf_(" maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations); - break; - } - } - - /*clean-up*/ - if(conserve_loads) delete loads; - VecFree(&uf); - VecFree(&ug); - VecFree(&old_uf); - VecFree(&old_ug); - -} Index: /issm/trunk/src/c/solvers/solvers.h =================================================================== --- /issm/trunk/src/c/solvers/solvers.h (revision 5996) +++ /issm/trunk/src/c/solvers/solvers.h (revision 5997) @@ -14,6 +14,4 @@ void solver_thermal_nonlinear(FemModel* femmodel); void solver_diagnostic_nonlinear(FemModel* femmodel,bool conserve_loads); -void solver_diagnostic_nonlinear_kgg(FemModel* femmodel,bool conserve_loads); -void solver_diagnostic_nonlinear_kff(FemModel* femmodel,bool conserve_loads); void solver_linear(FemModel* femmodel); void solver_adjoint_linear(FemModel* femmodel);