Index: /issm/trunk/src/c/modules/Reduceloadfromgtofx/Reduceloadfromgtofx.cpp =================================================================== --- /issm/trunk/src/c/modules/Reduceloadfromgtofx/Reduceloadfromgtofx.cpp (revision 5786) +++ /issm/trunk/src/c/modules/Reduceloadfromgtofx/Reduceloadfromgtofx.cpp (revision 5787) @@ -22,6 +22,10 @@ PetscScalar a; - /*Display message*/ - int verbose; parameters->FindParam(&verbose,VerboseEnum); + /*Parameters: */ + int verbose; + bool kffpartition=false; + + parameters->FindParam(&verbose,VerboseEnum); + parameters->FindParam(&kffpartition,KffEnum); if (verbose) _printf_(" Reducing Load vector from gset to fset\n"); @@ -39,5 +43,5 @@ if(nodesets->GetFSize()){ - VecPartition(&pf, pn, nodesets->GetPV_F(),nodesets->GetFSize()); + VecPartition(&pf, pn, nodesets->GetPV_F(),nodesets->GetFSize(),kffpartition); /*pf = pf - Kfs * y_s;*/ Index: /issm/trunk/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.cpp =================================================================== --- /issm/trunk/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.cpp (revision 5786) +++ /issm/trunk/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.cpp (revision 5787) @@ -5,9 +5,15 @@ #include "./Reducevectorgtofx.h" - -void Reducevectorgtofx(Vec* puf, Vec ug, NodeSets* nodesets){ + +void Reducevectorgtofx(Vec* puf, Vec ug, NodeSets* nodesets,Parameters* parameters){ /*output: */ Vec uf=NULL; + + /*what type of partitioning: */ + bool kffpartitioning=false; + + /*find parameter: */ + parameters->FindParam(&kffpartitioning,KffEnum); if(nodesets){ @@ -16,5 +22,5 @@ if (nodesets->GetGSize() && nodesets->GetFSize()){ - VecPartition(&uf,ug,nodesets->GetPV_F(),nodesets->GetFSize()); + VecPartition(&uf,ug,nodesets->GetPV_F(),nodesets->GetFSize(),kffpartitioning); } Index: /issm/trunk/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.h =================================================================== --- /issm/trunk/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.h (revision 5786) +++ /issm/trunk/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.h (revision 5787) @@ -10,5 +10,5 @@ /* local prototypes: */ -void Reducevectorgtofx(Vec* puf, Vec ug, NodeSets* nodesets); +void Reducevectorgtofx(Vec* puf, Vec ug, NodeSets* nodesets, Parameters* parameters); #endif /* _REDUCEVECTORGTOFX_H */ Index: /issm/trunk/src/c/modules/Reducevectorgtosx/Reducevectorgtosx.cpp =================================================================== --- /issm/trunk/src/c/modules/Reducevectorgtosx/Reducevectorgtosx.cpp (revision 5786) +++ /issm/trunk/src/c/modules/Reducevectorgtosx/Reducevectorgtosx.cpp (revision 5787) @@ -6,12 +6,18 @@ #include "./Reducevectorgtosx.h" -void Reducevectorgtosx( Vec* pys, Vec yg, NodeSets* nodesets){ +void Reducevectorgtosx( Vec* pys, Vec yg, NodeSets* nodesets,Parameters* parameters){ /*output: */ Vec ys=NULL; + /*what type of partitioning: */ + bool kffpartitioning=false; + + /*find parameter: */ + parameters->FindParam(&kffpartitioning,KffEnum); + if(nodesets){ if (nodesets->GetGSize() && nodesets->GetSSize()){ - VecPartition(&ys,yg,nodesets->GetPV_S(),nodesets->GetSSize()); + VecPartition(&ys,yg,nodesets->GetPV_S(),nodesets->GetSSize(),kffpartitioning); } } Index: /issm/trunk/src/c/modules/Reducevectorgtosx/Reducevectorgtosx.h =================================================================== --- /issm/trunk/src/c/modules/Reducevectorgtosx/Reducevectorgtosx.h (revision 5786) +++ /issm/trunk/src/c/modules/Reducevectorgtosx/Reducevectorgtosx.h (revision 5787) @@ -10,5 +10,5 @@ /* local prototypes: */ -void Reducevectorgtosx( Vec* pys, Vec yg, NodeSets* nodesets); +void Reducevectorgtosx( Vec* pys, Vec yg, NodeSets* nodesets,Parameters* parameters); #endif /* _REDUCEVECTORGTOSX_H */ Index: /issm/trunk/src/c/solutions/ResetBoundaryConditions.cpp =================================================================== --- /issm/trunk/src/c/solutions/ResetBoundaryConditions.cpp (revision 5786) +++ /issm/trunk/src/c/solutions/ResetBoundaryConditions.cpp (revision 5787) @@ -27,5 +27,5 @@ //Reduce from g to s set - Reducevectorgtosx(&ys,yg,femmodel->m_nodesets[analysis_counter]); + Reducevectorgtosx(&ys,yg,femmodel->m_nodesets[analysis_counter],femmodel->parameters); /*Plug into femmodel->m_ys: */ Index: /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear_kff.cpp =================================================================== --- /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear_kff.cpp (revision 5786) +++ /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear_kff.cpp (revision 5787) @@ -45,5 +45,5 @@ /*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); + 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) @@ -59,5 +59,5 @@ Reduceloadx(pf, Kfs, femmodel->ys, femmodel->parameters); MatFree(&Kfs); - + Solverx(&uf, Kff, pf, old_uf, femmodel->parameters); Index: /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear_kgg.cpp =================================================================== --- /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear_kgg.cpp (revision 5786) +++ /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear_kgg.cpp (revision 5787) @@ -45,5 +45,5 @@ /*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); + 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) Index: /issm/trunk/src/c/toolkits/petsc/patches/VecPartition.cpp =================================================================== --- /issm/trunk/src/c/toolkits/petsc/patches/VecPartition.cpp (revision 5786) +++ /issm/trunk/src/c/toolkits/petsc/patches/VecPartition.cpp (revision 5787) @@ -16,5 +16,5 @@ #include "../../../shared/shared.h" -void VecPartition(Vec* poutvector,Vec vectorA, double* row_partition_vector, int row_partition_vector_size){ +void VecPartition(Vec* poutvector,Vec vectorA, double* row_partition_vector, int row_partition_vector_size,bool kffpartition){ @@ -84,14 +84,10 @@ /*Ok, each node now have count values corresponding to the node_rows extracted values from vectorA. * From count and values, create a new vector*/ - - VecCreate(PETSC_COMM_WORLD,&outvector); ; + VecCreate(PETSC_COMM_WORLD,&outvector); VecSetSizes(outvector,count,PETSC_DECIDE); VecSetFromOptions(outvector); + VecGetOwnershipRange(outvector,&lower_row,&upper_row); - VecGetOwnershipRange(outvector,&lower_row,&upper_row); - upper_row--; - range=upper_row-lower_row+1; - - /*Set values using the values vector. First build new node_rows index.*/ + /*build new node_rows index.*/ if (count){ for (i=0;imax_nonlinear_iterations), - displaystring(femmodel.parameters.Verbose,'%s%i%s',' maximum number of iterations ',max_nonlinear_iterations,' exceeded'); - end - end - - end - - %deal with loads: - if conserve_loads==false, - femmodel.loads=loads; + if ~kffpartitioning, + femmodel=solver_diagnostic_nonlinear_kgg(femmodel,conserve_loads); + else + femmodel=solver_diagnostic_nonlinear_kff(femmodel,conserve_loads); end end Index: /issm/trunk/src/m/solvers/solver_diagnostic_nonlinear_kff.m =================================================================== --- /issm/trunk/src/m/solvers/solver_diagnostic_nonlinear_kff.m (revision 5787) +++ /issm/trunk/src/m/solvers/solver_diagnostic_nonlinear_kff.m (revision 5787) @@ -0,0 +1,73 @@ +function femmodel=solver_diagnostic_nonlinear_kff(femmodel,conserve_loads) +%SOLVER_DIAGNOSTIC_NONLINEAR - core solver of diagnostic run +% +% Usage: +% [femmodel]=solver_diagnostic_nonlinear_kff(femmodel,conserve_loads) + + %Recover parameters + min_mechanical_constraints=femmodel.parameters.MinMechanicalConstraints; + + %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,femmodel.parameters); + + %Update the solution to make sure that vx and vxold are similar + [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,ug); + + while(~converged), + + %save pointer to old velocity + old_ug=ug; + old_uf=uf; + + [K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters); + + p_f = Reduceload( p_f, K_fs, femmodel.ys); + + uf=Solver(K_ff,p_f,old_uf,femmodel.parameters); + + ug= Mergesolutionfromftog( uf, femmodel.ys, femmodel.nodesets,femmodel.parameters); + + [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,ug); + [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,femmodel.parameters); + + %add convergence status into status + [femmodel.elements loads]=InputUpdateFromConstant(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,double(converged),ConvergedEnum); + + %rift convergence + if ~constraints_converged, + if converged, + if num_unstable_constraints <= min_mechanical_constraints, + converged=1; + 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 + + %deal with loads: + if conserve_loads==false, + femmodel.loads=loads; + end +end Index: /issm/trunk/src/m/solvers/solver_diagnostic_nonlinear_kgg.m =================================================================== --- /issm/trunk/src/m/solvers/solver_diagnostic_nonlinear_kgg.m (revision 5787) +++ /issm/trunk/src/m/solvers/solver_diagnostic_nonlinear_kgg.m (revision 5787) @@ -0,0 +1,74 @@ +function femmodel=solver_diagnostic_nonlinear_kgg(femmodel,conserve_loads) +%SOLVER_DIAGNOSTIC_NONLINEAR - core solver of diagnostic run +% +% Usage: +% [femmodel]=solver_diagnostic_nonlinear_kgg(femmodel,conserve_loads) + + %Recover parameters + min_mechanical_constraints=femmodel.parameters.MinMechanicalConstraints; + + %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,femmodel.parameters); + + %Update the solution to make sure that vx and vxold are similar + [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,ug); + + while(~converged), + + %save pointer to old velocity + old_ug=ug; + old_uf=uf; + + [K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters); + + [K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.nodesets,femmodel.parameters); + p_f = Reduceloadfromgtof( p_g, K_fs, femmodel.ys, femmodel.nodesets,femmodel.parameters); + + uf=Solver(K_ff,p_f,old_uf,femmodel.parameters); + + ug= Mergesolutionfromftog( uf, femmodel.ys, femmodel.nodesets,femmodel.parameters); + + [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,ug); + [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,femmodel.parameters); + + %add convergence status into status + [femmodel.elements loads]=InputUpdateFromConstant(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,double(converged),ConvergedEnum); + + %rift convergence + if ~constraints_converged, + if converged, + if num_unstable_constraints <= min_mechanical_constraints, + converged=1; + 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 + + %deal with loads: + if conserve_loads==false, + femmodel.loads=loads; + end +end Index: /issm/trunk/src/mex/Makefile.am =================================================================== --- /issm/trunk/src/mex/Makefile.am (revision 5786) +++ /issm/trunk/src/mex/Makefile.am (revision 5787) @@ -49,4 +49,5 @@ ProcessParams\ Qmu\ + Reduceload\ Reduceloadfromgtof\ Reducematrixfromgtof\ @@ -231,4 +232,7 @@ Qmu/Qmu.h +Reduceload_SOURCES = Reduceload/Reduceload.cpp\ + Reduceload/Reduceload.h + Reducematrixfromgtof_SOURCES = Reducematrixfromgtof/Reducematrixfromgtof.cpp\ Reducematrixfromgtof/Reducematrixfromgtof.h Index: /issm/trunk/src/mex/Reduceload/Reduceload.cpp =================================================================== --- /issm/trunk/src/mex/Reduceload/Reduceload.cpp (revision 5787) +++ /issm/trunk/src/mex/Reduceload/Reduceload.cpp (revision 5787) @@ -0,0 +1,52 @@ +/*\file Reduceload.c + *\brief: reduce load from g set to f set + */ + +#include "./Reduceload.h" + +void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]){ + + /*input datasets: */ + Vec pf = NULL; + Mat Kfs = NULL; + Vec ys = NULL; + bool flag_ys0=false; + + /*Boot module: */ + MODULEBOOT(); + + /*checks on arguments on the matlab side: */ + if((nlhs!=NLHS) || (nrhs!=3 && nrhs!=4)){ + ReduceloadUsage(); + ISSMERROR(" usage. See above"); + } + + /*Input datasets: */ + FetchData(&pf,PF); + FetchData(&Kfs,KFS); + FetchData(&ys,YS); + if(nrhs==4)FetchData(&flag_ys0,YSFLAG); + else flag_ys0=false; + + /*!Reduce load from g to f size: */ + Reduceloadx(pf, Kfs, ys, flag_ys0); + + /*write output datasets: */ + WriteData(PFOUT,pf); + + /*Free ressources: */ + VecFree(&pf); + MatFree(&Kfs); + VecFree(&ys); + + /*end module: */ + MODULEEND(); +} + +void ReduceloadUsage(void) +{ + _printf_("\n"); + _printf_(" usage: [pf] = %s(pf,Kfs,ys);\n",__FUNCT__); + _printf_(" [pf] = %s(pf,Kfs,ys,ys0_flag);\n",__FUNCT__); + _printf_("\n"); +} Index: /issm/trunk/src/mex/Reduceload/Reduceload.h =================================================================== --- /issm/trunk/src/mex/Reduceload/Reduceload.h (revision 5787) +++ /issm/trunk/src/mex/Reduceload/Reduceload.h (revision 5787) @@ -0,0 +1,35 @@ + +/* + Reduceload.h +*/ + + +#ifndef _REDUCELOAD_H +#define _REDUCELOAD_H + +/* local prototypes: */ +void ReduceloadUsage(void); + +#include "../../c/modules/modules.h" +#include "../../c/Container/Container.h" +#include "../../c/shared/shared.h" + +#undef __FUNCT__ +#define __FUNCT__ "Reduceload" + +/* serial input macros: */ +#define PF (mxArray*)prhs[0] +#define KFS (mxArray*)prhs[1] +#define YS (mxArray*)prhs[2] +#define YSFLAG (mxArray*)prhs[3] + +/* serial output macros: */ +#define PFOUT (mxArray**)&plhs[0] + +/* serial arg counts: */ +#undef NLHS +#define NLHS 1 +#undef NRHS +#define NRHS 4 + +#endif /* _REDUCELOAD_H */ Index: /issm/trunk/src/mex/Reducevectorgtof/Reducevectorgtof.cpp =================================================================== --- /issm/trunk/src/mex/Reducevectorgtof/Reducevectorgtof.cpp (revision 5786) +++ /issm/trunk/src/mex/Reducevectorgtof/Reducevectorgtof.cpp (revision 5787) @@ -10,4 +10,5 @@ Vec ug=NULL; NodeSets* nodesets=NULL; + Parameters* parameters=NULL; /* output datasets: */ @@ -23,7 +24,8 @@ FetchData(&ug,UG); FetchNodeSets(&nodesets,NODESETS); + FetchParams(¶meters,PARAMETERS); /*!Reduce vector: */ - Reducevectorgtofx(&uf,ug,nodesets); + Reducevectorgtofx(&uf,ug,nodesets,parameters); /*write output datasets: */ @@ -32,4 +34,5 @@ /*Free ressources: */ delete nodesets; + delete parameters; VecFree(&ug); VecFree(&uf); @@ -42,5 +45,5 @@ { _printf_("\n"); - _printf_(" usage: uf = %s(ug,nodesets);\n",__FUNCT__); + _printf_(" usage: uf = %s(ug,nodesets,parameters);\n",__FUNCT__); _printf_("\n"); } Index: /issm/trunk/src/mex/Reducevectorgtof/Reducevectorgtof.h =================================================================== --- /issm/trunk/src/mex/Reducevectorgtof/Reducevectorgtof.h (revision 5786) +++ /issm/trunk/src/mex/Reducevectorgtof/Reducevectorgtof.h (revision 5787) @@ -21,4 +21,5 @@ #define UG (mxArray*)prhs[0] #define NODESETS (mxArray*)prhs[1] +#define PARAMETERS (mxArray*)prhs[2] /* serial output macros: */ @@ -29,5 +30,5 @@ #define NLHS 1 #undef NRHS -#define NRHS 2 +#define NRHS 3 Index: /issm/trunk/src/mex/Reducevectorgtos/Reducevectorgtos.cpp =================================================================== --- /issm/trunk/src/mex/Reducevectorgtos/Reducevectorgtos.cpp (revision 5786) +++ /issm/trunk/src/mex/Reducevectorgtos/Reducevectorgtos.cpp (revision 5787) @@ -10,4 +10,5 @@ Vec yg=NULL; NodeSets* nodesets=NULL; + Parameters* parameters=NULL; /* output datasets: */ @@ -23,7 +24,8 @@ FetchData(&yg,YG); FetchNodeSets(&nodesets,NODESETS); + FetchParams(¶meters,PARAMETERS); /*!Reduce vector: */ - Reducevectorgtosx(&ys,yg,nodesets); + Reducevectorgtosx(&ys,yg,nodesets,parameters); /*write output datasets: */ @@ -32,4 +34,5 @@ /*Free ressources: */ delete nodesets; + delete parameters; VecFree(&yg); VecFree(&ys); @@ -42,5 +45,5 @@ { _printf_("\n"); - _printf_(" usage: ys = %s(yg,nodesets);\n",__FUNCT__); + _printf_(" usage: ys = %s(yg,nodesets,parameters);\n",__FUNCT__); _printf_("\n"); } Index: /issm/trunk/src/mex/Reducevectorgtos/Reducevectorgtos.h =================================================================== --- /issm/trunk/src/mex/Reducevectorgtos/Reducevectorgtos.h (revision 5786) +++ /issm/trunk/src/mex/Reducevectorgtos/Reducevectorgtos.h (revision 5787) @@ -19,4 +19,5 @@ #define YG (mxArray*)prhs[0] #define NODESETS (mxArray*)prhs[1] +#define PARAMETERS (mxArray*)prhs[2] /* serial output macros: */ @@ -27,5 +28,5 @@ #define NLHS 1 #undef NRHS -#define NRHS 2 +#define NRHS 3 #endif /* _REDUCEVECTORGTOS_H */ Index: /issm/trunk/test/NightlyRun/test101.m =================================================================== --- /issm/trunk/test/NightlyRun/test101.m (revision 5786) +++ /issm/trunk/test/NightlyRun/test101.m (revision 5787) @@ -1,8 +1,12 @@ -md=mesh(model,'../Exp/Square.exp',150000); +md=mesh(model,'../Exp/Square.exp',50000); md=geography(md,'all',''); md=parameterize(md,'../Par/SquareShelfConstrained.par'); md=setelementstype(md,'macayeal','all'); md.cluster='none'; +md.kff=1; +tic md=solve(md,'analysis_type',DiagnosticSolutionEnum); +toc +error; %Fields and tolerances to track changes Index: /issm/trunk/test/NightlyRun/test102.m =================================================================== --- /issm/trunk/test/NightlyRun/test102.m (revision 5786) +++ /issm/trunk/test/NightlyRun/test102.m (revision 5787) @@ -4,4 +4,8 @@ md=setelementstype(md,'macayeal','all'); md.cluster=oshostname; +md.np=2; +md.verbose=1; +md.mem_debug=0; +md.kff=1; md=solve(md,'analysis_type',DiagnosticSolutionEnum);