/*!\file: solver_nonlinear.cpp * \brief: core of a non-linear solution, using fixed-point method */ #include "../toolkits/toolkits.h" #include "../objects/objects.h" #include "../io/io.h" #include "../EnumDefinitions/EnumDefinitions.h" #include "../modules/modules.h" #include "../solutions/solutions.h" #include "./solvers.h" void solver_newton(FemModel* femmodel){ /*intermediary: */ Mat Kff = NULL, Kfs = NULL, Jff = NULL; Vec ug = NULL, old_ug = NULL; Vec uf = NULL, old_uf = NULL, duf = NULL; Vec pf = NULL, pJf = NULL; Vec df = NULL; Vec ys = NULL; bool converged; int num_unstable_constraints; int count; /*parameters:*/ int max_nonlinear_iterations; int configuration_type; /*Recover parameters: */ femmodel->parameters->FindParam(&max_nonlinear_iterations,DiagnosticMaxiterEnum); femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); UpdateConstraintsx(femmodel->nodes,femmodel->constraints,femmodel->parameters); count=1; converged=false; /*Start non-linear iteration using input velocity: */ GetSolutionFromInputsx(&ug,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); Reducevectorgtofx(&uf,ug,femmodel->nodes,femmodel->parameters); //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate) InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,converged,ConvergedEnum); 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(&Kff,&Kfs,&pf,&df,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); Reduceloadx(pf,Kfs,ys);MatFree(&Kfs); Solverx(&uf,Kff,pf,old_uf,df,femmodel->parameters); convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters); CreateJacobianMatrixx(&Jff,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters); InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug); SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); VecDuplicate(pf,&pJf); MatMultPatch(Kff,uf,pJf); VecAXPY(pJf,-1.,pf); Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters); MatFree(&Kff);VecFree(&pf); VecFree(&df); VecAXPY(uf,-1.,duf); Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);VecFree(&ys); InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug); /*Increase count: */ count++; if(converged==true)break; if(count>=max_nonlinear_iterations){ _printf_(true," maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations); break; } } _printf_(VerboseConvergence(),"\n total number of iterations: %i\n",count-1); /*clean-up*/ VecFree(&uf); VecFree(&ug); VecFree(&old_ug); VecFree(&old_uf); }