Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17552) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17553) @@ -2900,9 +2900,7 @@ IssmDouble *xyz_list = NULL; - /*FIXME this should change*/ - r = 1.; - /*Get problem dimension*/ element->FindParam(&meshtype,MeshTypeEnum); + element->FindParam(&r,AugmentedLagrangianREnum); switch(meshtype){ case Mesh2DverticalEnum: dim = 2; epssize = 3; break; @@ -3289,5 +3287,5 @@ ElementVector* petemp =new ElementVector(pe1,pe2,pe3); ElementVector* pe4=CreatePVectorFSViscousXTH(element); - ElementVector* pe = new ElementVector(petemp,pe4); + pe = new ElementVector(petemp,pe4); delete pe1; delete pe2; @@ -3403,7 +3401,4 @@ } - /*FIXME*/ - r = 1.; - /*Fetch number of nodes and dof for this finite element*/ int vnumnodes = element->NumberofNodesVelocity(); @@ -3427,4 +3422,5 @@ /*Retrieve all inputs and parameters*/ + element->FindParam(&r,AugmentedLagrangianREnum); element->GetVerticesCoordinates(&xyz_list); @@ -4260,4 +4256,5 @@ Gauss* gauss = NULL; + parameters->FindParam(&r,AugmentedLagrangianREnum); parameters->FindParam(&meshtype,MeshTypeEnum); switch(meshtype){ @@ -4267,7 +4264,4 @@ default: _error_("mesh "<Size();i++){ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 17552) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 17553) @@ -624,4 +624,5 @@ AbsoluteEnum, IncrementalEnum, + AugmentedLagrangianREnum, /*}}}*/ /*Grounding Line{{{*/ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 17552) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 17553) @@ -603,4 +603,5 @@ case AbsoluteEnum : return "Absolute"; case IncrementalEnum : return "Incremental"; + case AugmentedLagrangianREnum : return "AugmentedLagrangianR"; case AgressiveMigrationEnum : return "AgressiveMigration"; case NoneEnum : return "None"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 17552) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 17553) @@ -615,4 +615,5 @@ else if (strcmp(name,"Absolute")==0) return AbsoluteEnum; else if (strcmp(name,"Incremental")==0) return IncrementalEnum; + else if (strcmp(name,"AugmentedLagrangianR")==0) return AugmentedLagrangianREnum; else if (strcmp(name,"AgressiveMigration")==0) return AgressiveMigrationEnum; else if (strcmp(name,"None")==0) return NoneEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"FSSolver")==0) return FSSolverEnum; else if (strcmp(name,"Adjoint")==0) return AdjointEnum; - else if (strcmp(name,"Colinear")==0) return ColinearEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"ControlSteady")==0) return ControlSteadyEnum; + if (strcmp(name,"Colinear")==0) return ColinearEnum; + else if (strcmp(name,"ControlSteady")==0) return ControlSteadyEnum; else if (strcmp(name,"Fset")==0) return FsetEnum; else if (strcmp(name,"Gradient1")==0) return Gradient1Enum; Index: /issm/trunk-jpl/src/c/shared/Numerics/NewtonSolveDnorm.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Numerics/NewtonSolveDnorm.cpp (revision 17552) +++ /issm/trunk-jpl/src/c/shared/Numerics/NewtonSolveDnorm.cpp (revision 17553) @@ -1,4 +1,5 @@ #include #include "./types.h" +#include "../Exceptions/exceptions.h" int NewtonSolveDnorm(IssmDouble* pdnorm,IssmDouble c1,IssmDouble c2,IssmDouble c3,IssmDouble n){ @@ -12,6 +13,12 @@ */ + /*trivial solution*/ + if(c3==0.){ + *pdnorm =0.; + return 0; + } /*Intermediaries*/ + int counter = 0; IssmDouble s = (1.+n)/n; IssmDouble y2; @@ -24,15 +31,19 @@ /*Newton step*/ - y2 = y1 - (c1*pow(pow(10.,y1),s-1.) + c2*pow(10.,y1) - c3)/((s-1)*c1*log(10)*pow(pow(10.,y1),s-1) + c2*log(10.)*pow(10.,y1)); + y2 = y1 - (c1*pow(pow(10.,y1),s-1.) + c2*pow(10.,y1) - c3)/((s-1)*c1*log(10)*pow(pow(10.,y1),s-1.) + c2*log(10.)*pow(10.,y1)); - if( fabs(y2-y1)/fabs(y1)50) _error_("Could not converge"); } /*Assign output pointer*/ *pdnorm = pow(10.,y2); + return 0; } Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp (revision 17552) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp (revision 17553) @@ -30,20 +30,31 @@ analysis->InitializeXTH(femmodel->elements,femmodel->parameters); - /*Solve KU=F*/ - SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); - CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); - Reduceloadx(pf, Kfs, ys); delete Kfs; - Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); - delete Kff; delete pf; delete df; - Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys; + /*Convergence criterion*/ + int count = 0; + femmodel->parameters->SetParam(10.,AugmentedLagrangianREnum); - /*Update solution*/ - InputUpdateFromSolutionx(femmodel,ug); + while(true){ - /*Update d and tau accordingly*/ - analysis->InputUpdateFromSolutionFSXTH(femmodel->elements,femmodel->parameters); + /*Solve KU=F*/ + SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); + CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); + Reduceloadx(pf, Kfs, ys); delete Kfs; + Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); + delete Kff; delete pf; delete df; + Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys; + + /*Update solution*/ + printf("ug norm = %g\n",ug->Norm(NORM_TWO)); + InputUpdateFromSolutionx(femmodel,ug); + + /*Update d and tau accordingly*/ + analysis->InputUpdateFromSolutionFSXTH(femmodel->elements,femmodel->parameters); + + count++; + if(count>5) break; + } /*Check for convergence*/ - _error_("STOP"); + //_error_("STOP"); delete ug; Index: /issm/trunk-jpl/src/m/enum/AugmentedLagrangianREnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/AugmentedLagrangianREnum.m (revision 17553) +++ /issm/trunk-jpl/src/m/enum/AugmentedLagrangianREnum.m (revision 17553) @@ -0,0 +1,11 @@ +function macro=AugmentedLagrangianREnum() +%AUGMENTEDLAGRANGIANRENUM - Enum of AugmentedLagrangianR +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=AugmentedLagrangianREnum() + +macro=StringToEnum('AugmentedLagrangianR'); Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 17552) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 17553) @@ -595,4 +595,5 @@ def AbsoluteEnum(): return StringToEnum("Absolute")[0] def IncrementalEnum(): return StringToEnum("Incremental")[0] +def AugmentedLagrangianREnum(): return StringToEnum("AugmentedLagrangianR")[0] def AgressiveMigrationEnum(): return StringToEnum("AgressiveMigration")[0] def NoneEnum(): return StringToEnum("None")[0]