Index: /issm/trunk/src/c/modules/Reducematrixfromgtofx/Reducematrixfromgtofx.cpp =================================================================== --- /issm/trunk/src/c/modules/Reducematrixfromgtofx/Reducematrixfromgtofx.cpp (revision 3907) +++ /issm/trunk/src/c/modules/Reducematrixfromgtofx/Reducematrixfromgtofx.cpp (revision 3907) @@ -0,0 +1,37 @@ +/*!\file Reducematrixfromgtofx + * \brief reduce matrix from g set to f fset + */ + +#include "./Reducematrixfromgtofx.h" + +void Reducematrixfromgtofx( Mat* pKff, Mat* pKfs,Mat Kgg,Mat Gmn,NodeSets* nodesets){ + + /*output: */ + Mat Kff=NULL; + Mat Kfs=NULL; + + /*intermediary: */ + Mat Knn=NULL; + + /*Reduce to n set: */ + Reducematrixfromgton( &Knn, Kgg, Gmn, nodesets->GetPV_M(), nodesets->GetMSize(),nodesets->GetPV_N(),nodesets->GetNSize(), 2 ); + + //Reduce matrix from n-size to f-size + if(nodesets->GetSSize()){ + + MatPartition(&Kff, Knn, nodesets->GetPV_F(),nodesets->GetFSize(),nodesets->GetPV_F(),nodesets->GetFSize()); + MatPartition(&Kfs, Knn, nodesets->GetPV_F(),nodesets->GetFSize(),nodesets->GetPV_S(),nodesets->GetSSize()); + + } + else{ + MatDuplicate(Knn,MAT_COPY_VALUES,&Kff); + Kfs=NULL; + } + + /*Free ressources:*/ + MatFree(&Knn); + + /*Assign output pointers:*/ + *pKff=Kff; + *pKfs=Kfs; +} Index: /issm/trunk/src/c/modules/Reducematrixfromgtofx/Reducematrixfromgtofx.h =================================================================== --- /issm/trunk/src/c/modules/Reducematrixfromgtofx/Reducematrixfromgtofx.h (revision 3907) +++ /issm/trunk/src/c/modules/Reducematrixfromgtofx/Reducematrixfromgtofx.h (revision 3907) @@ -0,0 +1,16 @@ +/*!\file: Reducematrixfromgtofx.h + * \brief reduce petsc matrix from g set to fset: + */ + + +#ifndef _REDUCEMATRIXFROMGTOFX_H +#define _REDUCEMATRIXFROMGTOFX_H + +#include "../objects/objects.h" + +/* local prototypes: */ +void Reducematrixfromgtofx( Mat* pKff, Mat* pKfs,Mat Kgg,Mat Gmn,NodeSets* nodesets); +void Reducematrixfromgton(Mat* pKnn,Mat Kgg,Mat Gmn,double* pv_m,int msize, double* pv_n,int nsize,int flag); + +#endif /* _REDUCEMATRIXFROMGTOFX_H */ + Index: /issm/trunk/src/c/modules/Reducematrixfromgtofx/Reducematrixfromgton.cpp =================================================================== --- /issm/trunk/src/c/modules/Reducematrixfromgtofx/Reducematrixfromgton.cpp (revision 3907) +++ /issm/trunk/src/c/modules/Reducematrixfromgtofx/Reducematrixfromgton.cpp (revision 3907) @@ -0,0 +1,73 @@ +/*!\file: Reducematrixfromgton + * \brief + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../toolkits/toolkits.h" + +void Reducematrixfromgton(Mat* pKnn,Mat Kgg,Mat Gmn,double* pv_m,int msize, double* pv_n,int nsize,int flag){ + + + /*output: */ + Mat Knn=NULL; + + /*intermediary: */ + Mat Kmm=NULL; + Mat Kmn=NULL; + Mat Knm=NULL; + Mat KnmGmn=NULL; + Mat KmmGmn=NULL; + Mat tGmn=NULL; + Mat tGmnKmn=NULL; + PetscScalar a=1; + + if (msize){ + + /*Partition K_gg*/ + MatPartition(&Kmm, Kgg, pv_m,msize,pv_m,msize); + MatPartition(&Kmn, Kgg, pv_m,msize,pv_n,nsize); + MatPartition(&Knm, Kgg, pv_n,nsize,pv_m,msize); + MatPartition(&Knn, Kgg, pv_n,nsize,pv_n,nsize); + + /*Reduce K_gg to K_nn*/ + //K_nn = K_nn + K_nm * G_mn; + MatMatMult(Knm,Gmn,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&KnmGmn); + MatAXPY(Knn,a,KnmGmn,DIFFERENT_NONZERO_PATTERN); + + if (flag!=1){ + //K_mn = K_mn + K_mm * G_mn; + MatMatMult(Kmm,Gmn,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&KmmGmn); + MatAXPY(Kmn,a,KmmGmn,DIFFERENT_NONZERO_PATTERN); + + //K_nn = K_nn + G_mn' * K_mn; + #if _PETSC_VERSION_ == 2 + MatTranspose(Gmn,&tGmn); + #else + MatTranspose(Gmn,MAT_INITIAL_MATRIX,&tGmn); + #endif + MatMatMult(tGmn,Kmn,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&tGmnKmn); + MatAXPY(Knn,a,tGmnKmn,DIFFERENT_NONZERO_PATTERN); + } + } + else{ + /*msize=0 just return input*/ + MatDuplicate(Kgg,MAT_COPY_VALUES,&Knn); + } + + /*Free ressources:*/ + MatFree(&Kmm); + MatFree(&Kmn); + MatFree(&Knm); + MatFree(&KnmGmn); + MatFree(&KmmGmn); + MatFree(&tGmn); + MatFree(&tGmnKmn); + + /*Assign output pointers:*/ + *pKnn=Knn; +} Index: /issm/trunk/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.cpp =================================================================== --- /issm/trunk/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.cpp (revision 3907) +++ /issm/trunk/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.cpp (revision 3907) @@ -0,0 +1,33 @@ +/*!\file Reducevectorgtofx + * \brief reduce petsc vector from g set to s set (free dofs), using the nodeset partitioning + * vectors. + */ + +#include "./Reducevectorgtofx.h" + +void Reducevectorgtofx(Vec* puf, Vec ug, NodeSets* nodesets){ + + /*output: */ + Vec uf=NULL; + + /*intermediary: */ + Vec un=NULL; + + if(nodesets){ + + + if (nodesets->GetNSize() && nodesets->GetFSize()){ + + VecPartition(&un,ug,nodesets->GetPV_N(),nodesets->GetNSize()); + + VecPartition(&uf,un,nodesets->GetPV_F(),nodesets->GetFSize()); + + } + + /*Free ressources:*/ + VecFree(&un); + } + + /*Assign output pointers:*/ + *puf=uf; +} Index: /issm/trunk/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.h =================================================================== --- /issm/trunk/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.h (revision 3907) +++ /issm/trunk/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.h (revision 3907) @@ -0,0 +1,15 @@ +/*!\file: Reducevectorgtofx.h + * \brief reduce petsc vector from g set to f set (free dofs), using the nodeset partitioning + * vectors. + */ + +#ifndef _REDUCEVECTORGTOFX_H +#define _REDUCEVECTORGTOFX_H + +#include "../objects/objects.h" + +/* local prototypes: */ +void Reducevectorgtofx(Vec* puf, Vec ug, NodeSets* nodesets); + +#endif /* _REDUCEVECTORGTOFX_H */ + Index: /issm/trunk/src/c/modules/Reducevectorgtosx/Reducevectorgtosx.cpp =================================================================== --- /issm/trunk/src/c/modules/Reducevectorgtosx/Reducevectorgtosx.cpp (revision 3907) +++ /issm/trunk/src/c/modules/Reducevectorgtosx/Reducevectorgtosx.cpp (revision 3907) @@ -0,0 +1,43 @@ +/*!\file Reducevectorgtosx + * \brief reduce petsc vector from g set to s set (single point constraints), using the nodeset partitioning + * vectors. + */ + +#include "./Reducevectorgtosx.h" + +void Reducevectorgtosx( Vec* pys, Vec* pys0, Vec yg, NodeSets* nodesets){ + + /*output: */ + Vec ys=NULL; + Vec ys0=NULL; + + /*intermediary: */ + Vec yn=NULL; + + if(nodesets){ + + + if (nodesets->GetNSize() && nodesets->GetSSize()){ + + VecPartition(&yn,yg,nodesets->GetPV_N(),nodesets->GetNSize()); + + VecPartition(&ys,yn,nodesets->GetPV_S(),nodesets->GetSSize()); + + } + + /*Create ys0, full of 0: */ + if(ys){ + VecDuplicate(ys,&ys0); + VecSet(ys0,0.0); + VecAssemblyBegin(ys0); + VecAssemblyEnd(ys0); + } + + /*Free ressources:*/ + VecFree(&yn); + } + + /*Assign output pointers:*/ + *pys=ys; + *pys0=ys0; +} Index: /issm/trunk/src/c/modules/Reducevectorgtosx/Reducevectorgtosx.h =================================================================== --- /issm/trunk/src/c/modules/Reducevectorgtosx/Reducevectorgtosx.h (revision 3907) +++ /issm/trunk/src/c/modules/Reducevectorgtosx/Reducevectorgtosx.h (revision 3907) @@ -0,0 +1,15 @@ +/*!\file: Reducevectorgtosx.h + * \brief reduce petsc vector from g set to s set (single point constraints), using the nodeset partitioning + * vectors. + */ + +#ifndef _REDUCEVECTORGTOSX_H +#define _REDUCEVECTORGTOSX_H + +#include "../objects/objects.h" + +/* local prototypes: */ +void Reducevectorgtosx( Vec* pys, Vec* pys0,Vec yg, NodeSets* nodesets); + +#endif /* _REDUCEVECTORGTOSX_H */ + Index: /issm/trunk/src/c/modules/Solverx/Solverx.cpp =================================================================== --- /issm/trunk/src/c/modules/Solverx/Solverx.cpp (revision 3907) +++ /issm/trunk/src/c/modules/Solverx/Solverx.cpp (revision 3907) @@ -0,0 +1,106 @@ +/*!\file Solverx + * \brief solver + */ + +#include "./Solverx.h" + +#include "../shared/shared.h" +#include "../include/include.h" + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + + +void Solverx( Vec* puf, Mat Kff, Vec pf, Vec uf0, char* solver_string){ + + /*output: */ + Vec uf=NULL; + + /*intermediary: */ + int local_m,local_n; + + /*Solver*/ + KSP ksp=NULL; + PC pc=NULL; + int iteration_number; + PetscTruth flag; + int solver_type; + + /*If initial guess for solution exists, use it to create uf, otherwise, + * duplicate right hand side so that solution vector has same structure*/ + if(uf0){ + VecDuplicate(uf0,&uf); VecCopy(uf0,uf); + } + else{ + MatGetLocalSize(Kff,&local_m,&local_n);uf=NewVecFromLocalSize(local_n); + } + + /*Before preparing the solver, add options to the options database*/ + PetscOptionsInsertMultipleString(solver_string); + + /*Process solver_string to see if we are not using special types of external solvers: */ + PetscOptionsDetermineSolverType(&solver_type,solver_string); + + #if _PETSC_VERSION_ == 2 + if (solver_type==MUMPSPACKAGE_LU){ + /*Convert Kff to MATTAIJMUMPS: */ + MatConvert(Kff,MATAIJMUMPS,MAT_REUSE_MATRIX,&Kff); + } + if (solver_type==MUMPSPACKAGE_CHOL){ + /*Convert Kff to MATTSBAIJMUMPS: */ + MatConvert(Kff,MATSBAIJMUMPS,MAT_REUSE_MATRIX,&Kff); + } + if (solver_type==SPOOLESPACKAGE_LU){ + /*Convert Kff to MATTSBAIJMUMPS: */ + MatConvert(Kff,MATAIJSPOOLES,MAT_REUSE_MATRIX,&Kff); + } + if (solver_type==SPOOLESPACKAGE_CHOL){ + /*Convert Kff to MATTSBAIJMUMPS: */ + MatConvert(Kff,MATSBAIJSPOOLES,MAT_REUSE_MATRIX,&Kff); + } + if (solver_type==SUPERLUDISTPACKAGE){ + /*Convert Kff to MATTSBAIJMUMPS: */ + MatConvert(Kff,MATSUPERLU_DIST,MAT_REUSE_MATRIX,&Kff); + } + #endif + + /*Prepare solver*/ + KSPCreate(MPI_COMM_WORLD,&ksp); + KSPSetOperators(ksp,Kff,Kff,DIFFERENT_NONZERO_PATTERN); + KSPSetFromOptions(ksp); + + #if _PETSC_VERSION_ == 3 + /*specific solver?: */ + KSPGetPC(ksp,&pc); + if (solver_type==MUMPSPACKAGE_LU){ + PCFactorSetMatSolverPackage(pc,MAT_SOLVER_MUMPS); + } + #endif + + + /*If initial guess for solution, use it, except if we are using the MUMPS direct solver, where any initial + * guess will crash Petsc: */ + if (uf0){ + if( (solver_type!=MUMPSPACKAGE_LU) && (solver_type!=MUMPSPACKAGE_CHOL) && (solver_type!=SPOOLESPACKAGE_LU)&& (solver_type!=SPOOLESPACKAGE_CHOL)&& (solver_type!=SUPERLUDISTPACKAGE)){ + KSPSetInitialGuessNonzero(ksp,PETSC_TRUE); + } + } + + KSPSolve(ksp,pf,uf); + + /*Check convergence*/ + KSPGetIterationNumber(ksp,&iteration_number); + if (iteration_number<0){ + ISSMERROR("%s%i"," Solver diverged at iteration number: ",-iteration_number); + } + + /*Free ressources:*/ + KSPFree(&ksp); + + + /*Assign output pointers:*/ + *puf=uf; +} Index: /issm/trunk/src/c/modules/Solverx/Solverx.h =================================================================== --- /issm/trunk/src/c/modules/Solverx/Solverx.h (revision 3907) +++ /issm/trunk/src/c/modules/Solverx/Solverx.h (revision 3907) @@ -0,0 +1,15 @@ +/*!\file: Solverx.h + * \brief solver + */ + +#ifndef _SOLVERX_H +#define _SOLVERX_H + +#include "../objects/objects.h" + + +/* local prototypes: */ +void Solverx( Vec* puf, Mat Kff, Vec pf, Vec uf0, char* solver_string); + +#endif /* _SOLVERX_H */ +