Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 11678) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 11679) @@ -344,4 +344,5 @@ ./modules/ResetCoordinateSystemx/ResetCoordinateSystemx.cpp\ ./modules/Solverx/Solverx.cpp\ + ./modules/Solverx/SolverxPetsc.cpp\ ./modules/Solverx/Solverx.h\ ./modules/Solverx/DofTypesToIndexSet.cpp\ Index: /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp (revision 11679) @@ -10,5 +10,5 @@ #include "../../EnumDefinitions/EnumDefinitions.h" -void CreateJacobianMatrixx(Mat* pJff,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,double kmax){ +void CreateJacobianMatrixx(Matrix** pJff,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,double kmax){ int i,connectivity; @@ -17,5 +17,5 @@ Element *element = NULL; Load *load = NULL; - Mat Jff = NULL; + Matrix* Jff = NULL; /*Checks*/ @@ -29,5 +29,5 @@ /*Initialize Jacobian Matrix*/ - Jff=NewMat(fsize,fsize,connectivity,numberofdofspernode); + Jff=new Matrix(fsize,fsize,connectivity,numberofdofspernode); /*Create and assemble matrix*/ @@ -41,6 +41,5 @@ if(load->InAnalysis(configuration_type)) load->PenaltyCreateJacobianMatrix(Jff,kmax); } - MatAssemblyBegin(Jff,MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(Jff,MAT_FINAL_ASSEMBLY); + Jff->Assemble(); /*Assign output pointer*/ Index: /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h (revision 11678) +++ /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h (revision 11679) @@ -10,5 +10,5 @@ /* local prototypes: */ -void CreateJacobianMatrixx(Mat* pJff,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,double kmax); +void CreateJacobianMatrixx(Matrix** pJff,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,double kmax); #endif /* _CREATEJACOBIANMATRIXX_H */ Index: /issm/trunk-jpl/src/c/modules/CreateNodalConstraintsx/CreateNodalConstraintsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/CreateNodalConstraintsx/CreateNodalConstraintsx.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/modules/CreateNodalConstraintsx/CreateNodalConstraintsx.cpp (revision 11679) @@ -10,5 +10,5 @@ #include "../../EnumDefinitions/EnumDefinitions.h" -void CreateNodalConstraintsx( Vec* pys, Nodes* nodes,int configuration_type){ +void CreateNodalConstraintsx( Vector** pys, Nodes* nodes,int configuration_type){ int i; @@ -18,5 +18,5 @@ /*output: */ - Vec ys=NULL; + Vector* ys=NULL; /*figure out how many dofs we have: */ @@ -24,5 +24,5 @@ /*allocate:*/ - ys=NewVec(numberofdofs); + ys=new Vector(numberofdofs); /*go through all nodes, and for the ones corresponding to this configuration_type, fill the @@ -36,6 +36,5 @@ /*Assemble: */ - VecAssemblyBegin(ys); - VecAssemblyEnd(ys); + ys->Assemble(); /*Assign output pointers: */ Index: /issm/trunk-jpl/src/c/modules/CreateNodalConstraintsx/CreateNodalConstraintsx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/CreateNodalConstraintsx/CreateNodalConstraintsx.h (revision 11678) +++ /issm/trunk-jpl/src/c/modules/CreateNodalConstraintsx/CreateNodalConstraintsx.h (revision 11679) @@ -9,5 +9,5 @@ /* local prototypes: */ -void CreateNodalConstraintsx( Vec* pys, Nodes* nodes,int configuration_type); +void CreateNodalConstraintsx( Vector** pys, Nodes* nodes,int configuration_type); #endif /* _CREATENODALCONSTRAINTSX_H */ Index: /issm/trunk-jpl/src/c/modules/GetSolutionFromInputsx/GetSolutionFromInputsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/GetSolutionFromInputsx/GetSolutionFromInputsx.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/modules/GetSolutionFromInputsx/GetSolutionFromInputsx.cpp (revision 11679) @@ -9,5 +9,5 @@ #include "../../EnumDefinitions/EnumDefinitions.h" -void GetSolutionFromInputsx( Vec* psolution, Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials, Parameters* parameters){ +void GetSolutionFromInputsx( Vector** psolution, Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials, Parameters* parameters){ /*intermediary: */ @@ -19,5 +19,5 @@ /*output: */ - Vec solution=NULL; + Vector* solution=NULL; /*retrive parameters: */ @@ -29,5 +29,5 @@ /*Initialize solution: */ - solution=NewVec(gsize); + solution=new Vector(gsize); /*Go through elements and plug solution: */ @@ -38,6 +38,5 @@ /*Assemble vector: */ - VecAssemblyBegin(solution); - VecAssemblyEnd(solution); + solution->Assemble(); /*Assign output pointers:*/ Index: /issm/trunk-jpl/src/c/modules/GetSolutionFromInputsx/GetSolutionFromInputsx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/GetSolutionFromInputsx/GetSolutionFromInputsx.h (revision 11678) +++ /issm/trunk-jpl/src/c/modules/GetSolutionFromInputsx/GetSolutionFromInputsx.h (revision 11679) @@ -10,5 +10,5 @@ /* local prototypes: */ -void GetSolutionFromInputsx( Vec* psolution, Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials, Parameters* parameters); +void GetSolutionFromInputsx( Vector** psolution, Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials, Parameters* parameters); #endif /* _GETSOLUTIONFROMINPUTSXX_H */ Index: /issm/trunk-jpl/src/c/modules/InputUpdateFromSolutionx/InputUpdateFromSolutionx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/InputUpdateFromSolutionx/InputUpdateFromSolutionx.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/modules/InputUpdateFromSolutionx/InputUpdateFromSolutionx.cpp (revision 11679) @@ -9,10 +9,10 @@ #include "../../EnumDefinitions/EnumDefinitions.h" -void InputUpdateFromSolutionx( Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,Vec solution){ +void InputUpdateFromSolutionx( Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,Vector* solution){ double* serial_solution=NULL; /*Serialize solution, so that elements can index into it on every CPU: */ - VecToMPISerial(&serial_solution,solution); + serial_solution=solution->ToMPISerial(); /*Call overloaded form of InputUpdateFromSolutionx: */ Index: /issm/trunk-jpl/src/c/modules/InputUpdateFromSolutionx/InputUpdateFromSolutionx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/InputUpdateFromSolutionx/InputUpdateFromSolutionx.h (revision 11678) +++ /issm/trunk-jpl/src/c/modules/InputUpdateFromSolutionx/InputUpdateFromSolutionx.h (revision 11679) @@ -10,9 +10,9 @@ /* local prototypes: */ -void InputUpdateFromSolutionx( Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials, Parameters* parameters,Vec solution); +void InputUpdateFromSolutionx( Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials, Parameters* parameters,Vector* solution); void InputUpdateFromSolutionx( Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,double* solution); //with timestep -void InputUpdateFromSolutionx( Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials, Parameters* parameters,Vec solution,int timestep); +void InputUpdateFromSolutionx( Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials, Parameters* parameters,Vector* solution,int timestep); void InputUpdateFromSolutionx( Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,double* solution, int timestep); Index: /issm/trunk-jpl/src/c/modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp (revision 11679) @@ -7,8 +7,8 @@ #include "./Mergesolutionfromftogx.h" -void Mergesolutionfromftogx( Vec* pug, Vec uf, Vec ys, Nodes* nodes, Parameters* parameters, bool flag_ys0){ +void Mergesolutionfromftogx( Vector** pug, Vector* uf, Vector* ys, Nodes* nodes, Parameters* parameters, bool flag_ys0){ /*output: */ - Vec ug=NULL; + Vector* ug=NULL; /*intermediary: */ @@ -29,10 +29,10 @@ if(ssize){ if(flag_ys0){ - VecSet(ys,0.0); + ys->Set(0.0); } } /*initialize ug: */ - ug=NewVec(gsize); + ug=new Vector(gsize); /*Merge f set back into g set: */ Index: /issm/trunk-jpl/src/c/modules/Mergesolutionfromftogx/Mergesolutionfromftogx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/Mergesolutionfromftogx/Mergesolutionfromftogx.h (revision 11678) +++ /issm/trunk-jpl/src/c/modules/Mergesolutionfromftogx/Mergesolutionfromftogx.h (revision 11679) @@ -9,5 +9,5 @@ /* local prototypes: */ -void Mergesolutionfromftogx( Vec* pug, Vec uf, Vec ys, Nodes* nodes, Parameters* parameters, bool flag_ys0=false); +void Mergesolutionfromftogx( Vector** pug, Vector* uf, Vector* ys, Nodes* nodes, Parameters* parameters, bool flag_ys0=false); #endif /* _MERGESOLUTIONFROMFTOGX_H */ Index: /issm/trunk-jpl/src/c/modules/Reduceloadx/Reduceloadx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Reduceloadx/Reduceloadx.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/modules/Reduceloadx/Reduceloadx.cpp (revision 11679) @@ -12,12 +12,11 @@ #include "../../io/io.h" -void Reduceloadx( Vec pf, Mat Kfs, Vec y_s,bool flag_ys0){ +void Reduceloadx( Vector* pf, Matrix* Kfs, Vector* y_s,bool flag_ys0){ /*intermediary*/ - Vec y_s0 = NULL; - Vec Kfsy_s = NULL; + Vector* y_s0 = NULL; + Vector* Kfsy_s = NULL; int Kfsm,Kfsn; int global_m,global_n; - PetscScalar a; bool fromlocalsize = true; int verbose; @@ -25,5 +24,5 @@ _printf_(VerboseModule()," Dirichlet lifting applied to load vector\n"); - MatGetSize(Kfs,&global_m,&global_n); + Kfs->GetSize(&global_m,&global_n); if(pf && global_m*global_n){ @@ -32,28 +31,26 @@ /*pf = pf - Kfs * y_s;*/ - MatGetLocalSize(Kfs,&Kfsm,&Kfsn); - Kfsy_s=NewVec(Kfsm,fromlocalsize); + Kfs->GetLocalSize(&Kfsm,&Kfsn); + Kfsy_s=new Vector(Kfsm,fromlocalsize); if (flag_ys0){ /*Create y_s0, full of 0: */ - VecDuplicate(y_s,&y_s0); - VecSet(y_s0,0.0); - VecAssemblyBegin(y_s0); - VecAssemblyEnd(y_s0); + y_s0=y_s->Duplicate(); + y_s0->Set(0.0); + y_s0->Assemble(); - MatMultPatch(Kfs,y_s0,Kfsy_s); + Kfs->MatMult(y_s0,Kfsy_s); } else{ - MatMultPatch(Kfs,y_s,Kfsy_s); + Kfs->MatMult(y_s,Kfsy_s); } - a=-1; - VecAXPY(pf,a,Kfsy_s); + pf->AXPY(Kfsy_s,-1); } /*Free ressources and return*/ - VecFree(&y_s0); - VecFree(&Kfsy_s); + delete y_s0; + delete Kfsy_s; } Index: /issm/trunk-jpl/src/c/modules/Reduceloadx/Reduceloadx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/Reduceloadx/Reduceloadx.h (revision 11678) +++ /issm/trunk-jpl/src/c/modules/Reduceloadx/Reduceloadx.h (revision 11679) @@ -9,5 +9,5 @@ /* local prototypes: */ -void Reduceloadx( Vec pf, Mat Kfs, Vec ys,bool flag_ys0=false); +void Reduceloadx( Vector* pf, Matrix* Kfs, Vector* ys,bool flag_ys0=false); #endif /* _REDUCELOADX_H */ Index: /issm/trunk-jpl/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.cpp (revision 11679) @@ -6,8 +6,8 @@ #include "./Reducevectorgtofx.h" -void Reducevectorgtofx(Vec* puf, Vec ug, Nodes* nodes,Parameters* parameters){ +void Reducevectorgtofx(Vector** puf, Vector* ug, Nodes* nodes,Parameters* parameters){ /*output: */ - Vec uf=NULL; + Vector* uf=NULL; /*variables: */ @@ -26,10 +26,10 @@ else{ /*allocate: */ - uf=NewVec(fsize); + uf=new Vector(fsize); if(nodes->NumberOfNodes(configuration_type)){ /*serialize ug, so nodes can index into it: */ - VecToMPISerial(&ug_serial,ug); + ug_serial=ug->ToMPISerial(); /*Go through all nodes, and ask them to retrieve values from ug, and plug them into uf: */ @@ -47,6 +47,5 @@ } /*Assemble vector: */ - VecAssemblyBegin(uf); - VecAssemblyEnd(uf); + uf->Assemble(); } Index: /issm/trunk-jpl/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.h (revision 11678) +++ /issm/trunk-jpl/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.h (revision 11679) @@ -10,5 +10,5 @@ /* local prototypes: */ -void Reducevectorgtofx(Vec* puf, Vec ug, Nodes* nodes,Parameters* parameters); +void Reducevectorgtofx(Vector** puf, Vector* ug, Nodes* nodes,Parameters* parameters); #endif /* _REDUCEVECTORGTOFX_H */ Index: /issm/trunk-jpl/src/c/modules/Reducevectorgtosx/Reducevectorgtosx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Reducevectorgtosx/Reducevectorgtosx.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/modules/Reducevectorgtosx/Reducevectorgtosx.cpp (revision 11679) @@ -6,8 +6,8 @@ #include "./Reducevectorgtosx.h" -void Reducevectorgtosx(Vec* pys, Vec yg, Nodes* nodes,Parameters* parameters){ +void Reducevectorgtosx(Vector** pys, Vector* yg, Nodes* nodes,Parameters* parameters){ /*output: */ - Vec ys=NULL; + Vector* ys=NULL; /*variables: */ @@ -26,10 +26,10 @@ else{ /*allocate: */ - ys=NewVec(ssize); + ys=new Vector(ssize); if(nodes->NumberOfNodes(configuration_type)){ /*serialize yg, so nodes can index into it: */ - VecToMPISerial(&yg_serial,yg); + yg_serial=yg->ToMPISerial(); /*Go throygh all nodes, and ask them to retrieve values from yg, and plyg them into ys: */ @@ -47,6 +47,5 @@ } /*Assemble vector: */ - VecAssemblyBegin(ys); - VecAssemblyEnd(ys); + ys->Assemble(); } Index: /issm/trunk-jpl/src/c/modules/Reducevectorgtosx/Reducevectorgtosx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/Reducevectorgtosx/Reducevectorgtosx.h (revision 11678) +++ /issm/trunk-jpl/src/c/modules/Reducevectorgtosx/Reducevectorgtosx.h (revision 11679) @@ -10,5 +10,5 @@ /* local prototypes: */ -void Reducevectorgtosx(Vec* pys, Vec yg, Nodes* nodes,Parameters* parameters); +void Reducevectorgtosx(Vector** pys, Vector* yg, Nodes* nodes,Parameters* parameters); #endif /* _REDUCEVECTORGTOSX_H */ Index: /issm/trunk-jpl/src/c/modules/Solverx/Solverx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Solverx/Solverx.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/modules/Solverx/Solverx.cpp (revision 11679) @@ -14,167 +14,21 @@ #endif -void Solverx(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters){ +void Solverx(Vector** puf, Matrix* Kff, Vector* pf, Vector* uf0,Vector* df, Parameters* parameters){ - /*Output: */ - Vec uf = NULL; + /*output: */ + Vector* uf=NULL; + uf=new Vector(); - /*Intermediary: */ - int local_m,local_n,global_m,global_n; - int analysis_type; + #ifdef _HAVE_PETSC_ + Vec uf0_vector=NULL; + if (uf0)uf0_vector=uf0->vector; - /*Solver */ - KSP ksp = NULL; - PC pc = NULL; - int iteration_number; - int solver_type; - bool fromlocalsize = true; - #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2) - PetscTruth flag,flg; + SolverxPetsc(&uf->vector, Kff->matrix, pf->vector, uf0_vector, df->vector, parameters); + VecGetSize(uf->vector,&uf->M); #else - PetscBool flag,flg; + _error_("not supported yet!"); #endif - /*Stokes: */ - IS isv=NULL; - IS isp=NULL; - - #if _PETSC_MAJOR_ >= 3 - char ksp_type[50]; - #endif - - - /*Display message*/ - _printf_(VerboseModule()," Solving\n"); - #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2) - if(VerboseSolver())PetscOptionsPrint(stdout); - #else - if(VerboseSolver())PetscOptionsView(PETSC_VIEWER_STDOUT_WORLD); - #endif - - /*First, check that f-set is not NULL, i.e. model is fully constrained: {{{*/ - _assert_(Kff); - MatGetSize(Kff,&global_m,&global_n); _assert_(global_m==global_m); - if(!global_n){ - *puf=NULL; return; - } - /*}}}*/ - /*Initial guess logic here: {{{1*/ - /*Now, check that we are not giving an initial guess to the solver, if we are running a direct solver: */ - #if _PETSC_MAJOR_ >= 3 - PetscOptionsGetString(PETSC_NULL,"-ksp_type",ksp_type,49,&flg); - if (strcmp(ksp_type,"preonly")==0)uf0=NULL; - #endif - - /*If initial guess for the solution exists, use it to create uf, otherwise, - * duplicate the right hand side so that the solution vector has the same structure*/ - if(uf0){ - VecDuplicate(uf0,&uf); VecCopy(uf0,uf); - } - else{ - MatGetLocalSize(Kff,&local_m,&local_n);uf=NewVec(local_n,fromlocalsize); - } - /*}}}*/ - /*Process petsc options to see if we are using special types of external solvers: {{{1*/ - PetscOptionsDetermineSolverType(&solver_type); - - /*In serial mode, the matrices have been loaded as MPIAIJ or AIJ matrices. - We need to convert them if we are going to run the solvers successfully: */ - #ifdef _SERIAL_ - #if _PETSC_MAJOR_ == 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); - } - if (solver_type==StokesSolverEnum){ - _error_("Petsc 2 does not support multi-physics solvers"); - } - #endif - #endif - /*}}}*/ - /*Check the solver is available: {{{1*/ - if(solver_type==MUMPSPACKAGE_LU || solver_type==MUMPSPACKAGE_CHOL){ - #if _PETSC_MAJOR_ >=3 - #ifndef _HAVE_MUMPS_ - _error_("requested MUMPS solver, which was not compiled into ISSM!\n"); - #endif - - #endif - } - /*}}}*/ - /*Prepare solver:{{{1*/ - KSPCreate(MPI_COMM_WORLD,&ksp); - KSPSetOperators(ksp,Kff,Kff,DIFFERENT_NONZERO_PATTERN); - KSPSetFromOptions(ksp); - - #if defined(_SERIAL_) && _PETSC_MAJOR_==3 - /*Specific solver?: */ - KSPGetPC(ksp,&pc); - if (solver_type==MUMPSPACKAGE_LU){ - #if _PETSC_MINOR_==1 - PCFactorSetMatSolverPackage(pc,MAT_SOLVER_MUMPS); - #else - PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS); - #endif - } - #endif - - #if defined(_PARALLEL_) && _PETSC_MAJOR_==3 - /*Stokes: */ - if (solver_type==StokesSolverEnum){ - /*Make indices out of doftypes: */ - if(!df)_error_("need doftypes for Stokes solver!\n"); - DofTypesToIndexSet(&isv,&isp,df,StokesSolverEnum); - - /*Set field splits: */ - KSPGetPC(ksp,&pc); - #if _PETSC_MINOR_==1 - PCFieldSplitSetIS(pc,isv); - PCFieldSplitSetIS(pc,isp); - #else - PCFieldSplitSetIS(pc,PETSC_NULL,isv); - PCFieldSplitSetIS(pc,PETSC_NULL,isp); - #endif - - } - #endif - - /*}}}*/ - /*If there is an initial guess for the solution, use it, except if we are using the MUMPS direct solver, where any initial guess will crash Petsc: {{{1*/ - 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); - } - } - /*}}}*/ - - if(VerboseSolver())KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD); - - /*Solve: */ - KSPSolve(ksp,pf,uf); - - /*Check convergence*/ - KSPGetIterationNumber(ksp,&iteration_number); - if (iteration_number<0) _error_("%s%i"," Solver diverged at iteration number: ",-iteration_number); - - /*Free resources:*/ - KSPFree(&ksp); - - /*Assign output pointers:*/ + /*Assign output pointers: */ *puf=uf; } Index: /issm/trunk-jpl/src/c/modules/Solverx/Solverx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/Solverx/Solverx.h (revision 11678) +++ /issm/trunk-jpl/src/c/modules/Solverx/Solverx.h (revision 11679) @@ -9,5 +9,6 @@ /* local prototypes: */ -void Solverx( Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df,Parameters* parameters); +void Solverx(Vector** puf, Matrix* Kff, Vector* pf, Vector* uf0,Vector* df, Parameters* parameters); +void SolverxPetsc(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters); void DofTypesToIndexSet(IS* pisv, IS* pisp, Vec df,int typeenum); #endif /* _SOLVERX_H */ Index: /issm/trunk-jpl/src/c/modules/Solverx/SolverxPetsc.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Solverx/SolverxPetsc.cpp (revision 11679) +++ /issm/trunk-jpl/src/c/modules/Solverx/SolverxPetsc.cpp (revision 11679) @@ -0,0 +1,180 @@ +/*!\file SolverxPetsc + * \brief Petsc implementation of solver + */ + +#include "./Solverx.h" +#include "../../shared/shared.h" +#include "../../include/include.h" +#include "../../io/io.h" + +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +void SolverxPetsc(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters){ + + /*Output: */ + Vec uf = NULL; + + /*Intermediary: */ + int local_m,local_n,global_m,global_n; + int analysis_type; + + /*Solver */ + KSP ksp = NULL; + PC pc = NULL; + int iteration_number; + int solver_type; + bool fromlocalsize = true; + #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2) + PetscTruth flag,flg; + #else + PetscBool flag,flg; + #endif + + /*Stokes: */ + IS isv=NULL; + IS isp=NULL; + + #if _PETSC_MAJOR_ >= 3 + char ksp_type[50]; + #endif + + + /*Display message*/ + _printf_(VerboseModule()," Solving\n"); + #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2) + if(VerboseSolver())PetscOptionsPrint(stdout); + #else + if(VerboseSolver())PetscOptionsView(PETSC_VIEWER_STDOUT_WORLD); + #endif + + /*First, check that f-set is not NULL, i.e. model is fully constrained: {{{*/ + _assert_(Kff); + MatGetSize(Kff,&global_m,&global_n); _assert_(global_m==global_m); + if(!global_n){ + *puf=NULL; return; + } + /*}}}*/ + /*Initial guess logic here: {{{1*/ + /*Now, check that we are not giving an initial guess to the solver, if we are running a direct solver: */ + #if _PETSC_MAJOR_ >= 3 + PetscOptionsGetString(PETSC_NULL,"-ksp_type",ksp_type,49,&flg); + if (strcmp(ksp_type,"preonly")==0)uf0=NULL; + #endif + + /*If initial guess for the solution exists, use it to create uf, otherwise, + * duplicate the right hand side so that the solution vector has the same structure*/ + if(uf0){ + VecDuplicate(uf0,&uf); VecCopy(uf0,uf); + } + else{ + MatGetLocalSize(Kff,&local_m,&local_n);uf=NewVec(local_n,fromlocalsize); + } + /*}}}*/ + /*Process petsc options to see if we are using special types of external solvers: {{{1*/ + PetscOptionsDetermineSolverType(&solver_type); + + /*In serial mode, the matrices have been loaded as MPIAIJ or AIJ matrices. + We need to convert them if we are going to run the solvers successfully: */ + #ifdef _SERIAL_ + #if _PETSC_MAJOR_ == 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); + } + if (solver_type==StokesSolverEnum){ + _error_("Petsc 2 does not support multi-physics solvers"); + } + #endif + #endif + /*}}}*/ + /*Check the solver is available: {{{1*/ + if(solver_type==MUMPSPACKAGE_LU || solver_type==MUMPSPACKAGE_CHOL){ + #if _PETSC_MAJOR_ >=3 + #ifndef _HAVE_MUMPS_ + _error_("requested MUMPS solver, which was not compiled into ISSM!\n"); + #endif + + #endif + } + /*}}}*/ + /*Prepare solver:{{{1*/ + KSPCreate(MPI_COMM_WORLD,&ksp); + KSPSetOperators(ksp,Kff,Kff,DIFFERENT_NONZERO_PATTERN); + KSPSetFromOptions(ksp); + + #if defined(_SERIAL_) && _PETSC_MAJOR_==3 + /*Specific solver?: */ + KSPGetPC(ksp,&pc); + if (solver_type==MUMPSPACKAGE_LU){ + #if _PETSC_MINOR_==1 + PCFactorSetMatSolverPackage(pc,MAT_SOLVER_MUMPS); + #else + PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS); + #endif + } + #endif + + #if defined(_PARALLEL_) && _PETSC_MAJOR_==3 + /*Stokes: */ + if (solver_type==StokesSolverEnum){ + /*Make indices out of doftypes: */ + if(!df)_error_("need doftypes for Stokes solver!\n"); + DofTypesToIndexSet(&isv,&isp,df,StokesSolverEnum); + + /*Set field splits: */ + KSPGetPC(ksp,&pc); + #if _PETSC_MINOR_==1 + PCFieldSplitSetIS(pc,isv); + PCFieldSplitSetIS(pc,isp); + #else + PCFieldSplitSetIS(pc,PETSC_NULL,isv); + PCFieldSplitSetIS(pc,PETSC_NULL,isp); + #endif + + } + #endif + + /*}}}*/ + /*If there is an initial guess for the solution, use it, except if we are using the MUMPS direct solver, where any initial guess will crash Petsc: {{{1*/ + 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); + } + } + /*}}}*/ + + if(VerboseSolver())KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD); + + /*Solve: */ + KSPSolve(ksp,pf,uf); + + /*Check convergence*/ + KSPGetIterationNumber(ksp,&iteration_number); + if (iteration_number<0) _error_("%s%i"," Solver diverged at iteration number: ",-iteration_number); + + /*Free resources:*/ + KSPFree(&ksp); + + /*Assign output pointers:*/ + *puf=uf; +} Index: /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp (revision 11679) @@ -10,5 +10,5 @@ #include "../../EnumDefinitions/EnumDefinitions.h" -void SystemMatricesx(Mat* pKff, Mat* pKfs, Vec* ppf, Vec* pdf, double* pkmax,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,bool kflag,bool pflag,bool penalty_kflag,bool penalty_pflag){ +void SystemMatricesx(Matrix** pKff, Matrix** pKfs, Vector** ppf, Vector** pdf, double* pkmax,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,bool kflag,bool pflag,bool penalty_kflag,bool penalty_pflag){ /*intermediary: */ @@ -21,8 +21,8 @@ /*output: */ - Mat Kff = NULL; - Mat Kfs = NULL; - Vec pf = NULL; - Vec df=NULL; + Matrix* Kff = NULL; + Matrix* Kfs = NULL; + Vector* pf = NULL; + Vector* df=NULL; double kmax = 0; @@ -49,7 +49,7 @@ if(kflag){ - Kff=NewMat(fsize,fsize,connectivity,numberofdofspernode); - Kfs=NewMat(fsize,ssize,connectivity,numberofdofspernode); - df=NewVec(fsize); + Kff=new Matrix(fsize,fsize,connectivity,numberofdofspernode); + Kfs=new Matrix(fsize,ssize,connectivity,numberofdofspernode); + df=new Vector(fsize); /*Fill stiffness matrix from elements: */ @@ -66,17 +66,8 @@ /*Assemble matrix and doftypes and compress matrix to save memory: */ - MatAssemblyBegin(Kff,MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(Kff,MAT_FINAL_ASSEMBLY); - #if _PETSC_MAJOR_ == 2 - MatCompress(Kff); - #endif + Kff->Assemble(); + Kfs->Assemble(); + df->Assemble(); - MatAssemblyBegin(Kfs,MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(Kfs,MAT_FINAL_ASSEMBLY); - #if _PETSC_MAJOR_ == 2 - MatCompress(Kfs); - #endif - VecAssemblyBegin(df); - VecAssemblyEnd(df); } @@ -84,5 +75,5 @@ if(pflag){ - pf=NewVec(fsize); + pf=new Vector(fsize); /*Fill right hand side vector, from elements: */ @@ -98,10 +89,9 @@ } - VecAssemblyBegin(pf); - VecAssemblyEnd(pf); + pf->Assemble(); } /*Now, figure out maximum value of K_gg, so that we can penalize it correctly: */ - MatNorm(Kff,NORM_INFINITY,&kmax); + kmax=Kff->Norm(NORM_INFINITY); /*Now, deal with penalties*/ @@ -115,15 +105,6 @@ /*Assemble matrix and compress matrix to save memory: */ - MatAssemblyBegin(Kff,MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(Kff,MAT_FINAL_ASSEMBLY); - #if _PETSC_MAJOR_ == 2 - MatCompress(Kff); - #endif - - MatAssemblyBegin(Kfs,MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(Kfs,MAT_FINAL_ASSEMBLY); - #if _PETSC_MAJOR_ == 2 - MatCompress(Kfs); - #endif + Kff->Assemble(); + Kfs->Assemble(); } @@ -137,17 +118,17 @@ } - VecAssemblyBegin(pf); - VecAssemblyEnd(pf); + pf->Assemble(); + pf->Assemble(); } /*Assign output pointers: */ if(pKff) *pKff=Kff; - else MatFree(&Kff); + else delete Kff; if(pKfs) *pKfs=Kfs; - else MatFree(&Kfs); + else delete Kfs; if(ppf) *ppf=pf; - else VecFree(&pf); + else delete pf; if(pdf) *pdf=df; - else VecFree(&df); + else delete df; if(pkmax) *pkmax=kmax; } Index: /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.h (revision 11678) +++ /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.h (revision 11679) @@ -10,5 +10,5 @@ /* local prototypes: */ -void SystemMatricesx(Mat* pKff, Mat* pKfs, Vec* ppf, Vec* pdf, double* pkmax,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters, +void SystemMatricesx(Matrix** pKff, Matrix** pKfs, Vector** ppf, Vector** pdf, double* pkmax,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters, bool kflag=true,bool pflag=true,bool penalty_kflag=true,bool penalty_pflag=true); Index: /issm/trunk-jpl/src/c/modules/UpdateDynamicConstraintsx/UpdateDynamicConstraintsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/UpdateDynamicConstraintsx/UpdateDynamicConstraintsx.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/modules/UpdateDynamicConstraintsx/UpdateDynamicConstraintsx.cpp (revision 11679) @@ -10,5 +10,5 @@ #include "../../EnumDefinitions/EnumDefinitions.h" -void UpdateDynamicConstraintsx(Constraints* constraints,Nodes* nodes,Parameters* parameters,Vec yg){ +void UpdateDynamicConstraintsx(Constraints* constraints,Nodes* nodes,Parameters* parameters,Vector* yg){ int configuration_type; @@ -19,5 +19,5 @@ /*serialize yg, so nodes can index into it: */ - VecToMPISerial(&yg_serial,yg); + yg_serial=yg->ToMPISerial(); for(int i=0;iSize();i++){ Index: /issm/trunk-jpl/src/c/modules/UpdateDynamicConstraintsx/UpdateDynamicConstraintsx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/UpdateDynamicConstraintsx/UpdateDynamicConstraintsx.h (revision 11678) +++ /issm/trunk-jpl/src/c/modules/UpdateDynamicConstraintsx/UpdateDynamicConstraintsx.h (revision 11679) @@ -9,5 +9,5 @@ #include "../../objects/objects.h" -void UpdateDynamicConstraintsx(Constraints* constraints,Nodes* nodes,Parameters* parameters,Vec yg); +void UpdateDynamicConstraintsx(Constraints* constraints,Nodes* nodes,Parameters* parameters,Vector* yg); #endif /* _UPDATESPCSX_H */ Index: /issm/trunk-jpl/src/c/modules/VecMergex/VecMergex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/VecMergex/VecMergex.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/modules/VecMergex/VecMergex.cpp (revision 11679) @@ -10,5 +10,5 @@ #include "../../EnumDefinitions/EnumDefinitions.h" -void VecMergex(Vec ug, Vec uf, Nodes* nodes, Parameters* parameters, int SetEnum){ +void VecMergex(Vector* ug, Vector* uf, Nodes* nodes, Parameters* parameters, int SetEnum){ /*variables: */ @@ -21,5 +21,6 @@ /*serialize uf: */ - VecToMPISerial(&uf_serial,uf); + uf_serial=uf->ToMPISerial(); + /*Do we have any nodes for this configuration? :*/ @@ -43,6 +44,4 @@ /*Assemble vector: */ - VecAssemblyBegin(ug); - VecAssemblyEnd(ug); - + ug->Assemble(); } Index: /issm/trunk-jpl/src/c/modules/VecMergex/VecMergex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/VecMergex/VecMergex.h (revision 11678) +++ /issm/trunk-jpl/src/c/modules/VecMergex/VecMergex.h (revision 11679) @@ -10,5 +10,5 @@ /* local prototypes: */ -void VecMergex(Vec ug, Vec uf, Nodes* nodes, Parameters* parameters, int SetEnum); +void VecMergex(Vector* ug, Vector* uf, Nodes* nodes, Parameters* parameters, int SetEnum); #endif /* _VECMERGEX_H */ Index: /issm/trunk-jpl/src/c/objects/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Element.h (revision 11678) +++ /issm/trunk-jpl/src/c/objects/Elements/Element.h (revision 11679) @@ -16,4 +16,6 @@ class Parameters; class Patch; +class Matrix; +class Vector; #include "../../toolkits/toolkits.h" @@ -28,8 +30,8 @@ virtual void Configure(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters)=0; virtual void SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters)=0; - virtual void CreateKMatrix(Mat Kff, Mat Kfs,Vec df)=0; - virtual void CreatePVector(Vec pf)=0; - virtual void CreateJacobianMatrix(Mat Jff)=0; - virtual void GetSolutionFromInputs(Vec solution)=0; + virtual void CreateKMatrix(Matrix* Kff, Matrix* Kfs,Vector* df)=0; + virtual void CreatePVector(Vector* pf)=0; + virtual void CreateJacobianMatrix(Matrix* Jff)=0; + virtual void GetSolutionFromInputs(Vector* solution)=0; virtual int GetNodeIndex(Node* node)=0; virtual int Sid()=0; Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp (revision 11679) @@ -567,5 +567,5 @@ /*}}}*/ /*FUNCTION Penta::CreateKMatrix {{{1*/ -void Penta::CreateKMatrix(Mat Kff, Mat Kfs,Vec df){ +void Penta::CreateKMatrix(Matrix* Kff, Matrix* Kfs,Vector* df){ /*retrieve parameters: */ @@ -671,5 +671,5 @@ /*}}}*/ /*FUNCTION Penta::CreatePVector {{{1*/ -void Penta::CreatePVector(Vec pf){ +void Penta::CreatePVector(Vector* pf){ /*retrive parameters: */ @@ -773,5 +773,5 @@ /*}}}*/ /*FUNCTION Penta::CreateJacobianMatrix{{{1*/ -void Penta::CreateJacobianMatrix(Mat Jff){ +void Penta::CreateJacobianMatrix(Matrix* Jff){ /*retrieve parameters: */ @@ -1091,5 +1091,5 @@ /*}}}*/ /*FUNCTION Penta::GetSolutionFromInputs{{{1*/ -void Penta::GetSolutionFromInputs(Vec solution){ +void Penta::GetSolutionFromInputs(Vector* solution){ int analysis_type; @@ -4180,5 +4180,5 @@ /*}}}*/ /*FUNCTION Penta::GetSolutionFromInputsThermal{{{1*/ -void Penta::GetSolutionFromInputsThermal(Vec solution){ +void Penta::GetSolutionFromInputsThermal(Vector* solution){ const int numdof=NDOF1*NUMVERTICES; @@ -4203,5 +4203,5 @@ /*Add value to global vector*/ - VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); + solution->SetValues(numdof,doflist,values,INSERT_VALUES); /*Free ressources:*/ @@ -4211,5 +4211,5 @@ /*}}}*/ /*FUNCTION Penta::GetSolutionFromInputsEnthalpy{{{1*/ -void Penta::GetSolutionFromInputsEnthalpy(Vec solution){ +void Penta::GetSolutionFromInputsEnthalpy(Vector* solution){ const int numdof=NDOF1*NUMVERTICES; @@ -4234,5 +4234,5 @@ /*Add value to global vector*/ - VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); + solution->SetValues(numdof,doflist,values,INSERT_VALUES); /*Free ressources:*/ @@ -7840,5 +7840,5 @@ /*}}}*/ /*FUNCTION Penta::GetSolutionFromInputsDiagnosticHoriz{{{1*/ -void Penta::GetSolutionFromInputsDiagnosticHoriz(Vec solution){ +void Penta::GetSolutionFromInputsDiagnosticHoriz(Vector* solution){ const int numdof=NDOF2*NUMVERTICES; @@ -7874,5 +7874,5 @@ /*Add value to global vector*/ - VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); + solution->SetValues(numdof,doflist,values,INSERT_VALUES); /*Free ressources:*/ @@ -7882,5 +7882,5 @@ /*}}}*/ /*FUNCTION Penta::GetSolutionFromInputsDiagnosticHutter{{{1*/ -void Penta::GetSolutionFromInputsDiagnosticHutter(Vec solution){ +void Penta::GetSolutionFromInputsDiagnosticHutter(Vector* solution){ const int numdof=NDOF2*NUMVERTICES; @@ -7910,5 +7910,5 @@ /*Add value to global vector*/ - VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); + solution->SetValues(numdof,doflist,values,INSERT_VALUES); /*Free ressources:*/ @@ -7918,5 +7918,5 @@ /*}}}*/ /*FUNCTION Penta::GetSolutionFromInputsDiagnosticVert{{{1*/ -void Penta::GetSolutionFromInputsDiagnosticVert(Vec solution){ +void Penta::GetSolutionFromInputsDiagnosticVert(Vector* solution){ const int numdof=NDOF1*NUMVERTICES; @@ -7943,5 +7943,5 @@ /*Add value to global vector*/ - VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); + solution->SetValues(numdof,doflist,values,INSERT_VALUES); /*Free ressources:*/ @@ -7951,5 +7951,5 @@ /*}}}*/ /*FUNCTION Penta::GetSolutionFromInputsDiagnosticStokes{{{1*/ -void Penta::GetSolutionFromInputsDiagnosticStokes(Vec solution){ +void Penta::GetSolutionFromInputsDiagnosticStokes(Vector* solution){ const int numdof=NDOF4*NUMVERTICES; @@ -7988,5 +7988,5 @@ /*Add value to global vector*/ - VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); + solution->SetValues(numdof,doflist,values,INSERT_VALUES); /*Free ressources:*/ Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Penta.h (revision 11678) +++ /issm/trunk-jpl/src/c/objects/Elements/Penta.h (revision 11679) @@ -86,10 +86,10 @@ void Configure(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters); void SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters); - void CreateKMatrix(Mat Kff, Mat Kfs,Vec df); - void CreatePVector(Vec pf); - void CreateJacobianMatrix(Mat Jff); + void CreateKMatrix(Matrix* Kff, Matrix* Kfs,Vector* df); + void CreatePVector(Vector* pf); + void CreateJacobianMatrix(Matrix* Jff); void DeleteResults(void); int GetNodeIndex(Node* node); - void GetSolutionFromInputs(Vec solution); + void GetSolutionFromInputs(Vector* solution); double GetZcoord(GaussPenta* gauss); void GetVectorFromInputs(Vec vector,int name_enum); @@ -183,5 +183,5 @@ void GetInputValue(double* pvalue,Node* node,int enumtype); void GetPhi(double* phi, double* epsilon, double viscosity); - void GetSolutionFromInputsEnthalpy(Vec solutiong); + void GetSolutionFromInputsEnthalpy(Vector* solutiong); double GetStabilizationParameter(double u, double v, double w, double diameter, double kappa); void GetStrainRate3dPattyn(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input); @@ -250,8 +250,8 @@ void InputUpdateFromSolutionDiagnosticVert( double* solutiong); void InputUpdateFromSolutionDiagnosticStokes( double* solutiong); - void GetSolutionFromInputsDiagnosticHoriz(Vec solutiong); - void GetSolutionFromInputsDiagnosticHutter(Vec solutiong); - void GetSolutionFromInputsDiagnosticStokes(Vec solutiong); - void GetSolutionFromInputsDiagnosticVert(Vec solutiong); + void GetSolutionFromInputsDiagnosticHoriz(Vector* solutiong); + void GetSolutionFromInputsDiagnosticHutter(Vector* solutiong); + void GetSolutionFromInputsDiagnosticStokes(Vector* solutiong); + void GetSolutionFromInputsDiagnosticVert(Vector* solutiong); ElementVector* CreatePVectorCouplingMacAyealStokes(void); ElementVector* CreatePVectorCouplingMacAyealStokesViscous(void); @@ -307,5 +307,5 @@ ElementVector* CreatePVectorThermalShelf(void); ElementVector* CreatePVectorThermalSheet(void); - void GetSolutionFromInputsThermal(Vec solutiong); + void GetSolutionFromInputsThermal(Vector* solutiong); void InputUpdateFromSolutionThermal( double* solutiong); void InputUpdateFromSolutionEnthalpy( double* solutiong); Index: /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp (revision 11679) @@ -319,5 +319,5 @@ /*}}}*/ /*FUNCTION Tria::CreateKMatrix {{{1*/ -void Tria::CreateKMatrix(Mat Kff, Mat Kfs,Vec df){ +void Tria::CreateKMatrix(Matrix* Kff, Matrix* Kfs,Vector* df){ /*retreive parameters: */ @@ -672,5 +672,5 @@ /*}}}*/ /*FUNCTION Tria::CreatePVector {{{1*/ -void Tria::CreatePVector(Vec pf){ +void Tria::CreatePVector(Vector* pf){ /*retrive parameters: */ @@ -892,5 +892,5 @@ /*}}}*/ /*FUNCTION Tria::CreateJacobianMatrix{{{1*/ -void Tria::CreateJacobianMatrix(Mat Jff){ +void Tria::CreateJacobianMatrix(Matrix* Jff){ /*retrieve parameters: */ @@ -1274,5 +1274,5 @@ /*}}}*/ /*FUNCTION Tria::GetSolutionFromInputs{{{1*/ -void Tria::GetSolutionFromInputs(Vec solution){ +void Tria::GetSolutionFromInputs(Vector* solution){ /*retrive parameters: */ @@ -3139,5 +3139,5 @@ /*}}}*/ /*FUNCTION Tria::GetSolutionFromInputsDiagnosticHoriz{{{1*/ -void Tria::GetSolutionFromInputsDiagnosticHoriz(Vec solution){ +void Tria::GetSolutionFromInputsDiagnosticHoriz(Vector* solution){ const int numdof=NDOF2*NUMVERTICES; @@ -3170,5 +3170,5 @@ } - VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); + solution->SetValues(numdof,doflist,values,INSERT_VALUES); /*Free ressources:*/ @@ -3178,5 +3178,5 @@ /*}}}*/ /*FUNCTION Tria::GetSolutionFromInputsDiagnosticHutter{{{1*/ -void Tria::GetSolutionFromInputsDiagnosticHutter(Vec solution){ +void Tria::GetSolutionFromInputsDiagnosticHutter(Vector* solution){ const int numdof=NDOF2*NUMVERTICES; @@ -3209,5 +3209,5 @@ } - VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); + solution->SetValues(numdof,doflist,values,INSERT_VALUES); /*Free ressources:*/ @@ -4881,4 +4881,5 @@ } /*}}}*/ + /*FUNCTION Tria::InputUpdateFromSolutionAdjointBalancethickness {{{1*/ void Tria::InputUpdateFromSolutionAdjointBalancethickness(double* solution){ @@ -5185,5 +5186,5 @@ /*}}}*/ /*FUNCTION Tria::GetSolutionFromInputsHydrology{{{1*/ -void Tria::GetSolutionFromInputsHydrology(Vec solution){ +void Tria::GetSolutionFromInputsHydrology(Vector* solution){ const int numdof=NDOF1*NUMVERTICES; @@ -5213,5 +5214,5 @@ } - VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); + solution->SetValues(numdof,doflist,values,INSERT_VALUES); /*Free ressources:*/ Index: /issm/trunk-jpl/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Tria.h (revision 11678) +++ /issm/trunk-jpl/src/c/objects/Elements/Tria.h (revision 11679) @@ -82,7 +82,7 @@ void Configure(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters); void SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters); - void CreateKMatrix(Mat Kff, Mat Kfs,Vec df); - void CreatePVector(Vec pf); - void CreateJacobianMatrix(Mat Jff); + void CreateKMatrix(Matrix* Kff, Matrix* Kfs,Vector* df); + void CreatePVector(Vector* pf); + void CreateJacobianMatrix(Matrix* Jff); int GetNodeIndex(Node* node); int Sid(); @@ -92,5 +92,5 @@ bool IsNodeOnShelfFromFlags(double* flags); bool IsOnWater(); - void GetSolutionFromInputs(Vec solution); + void GetSolutionFromInputs(Vector* solution); void GetVectorFromInputs(Vec vector, int name_enum); void GetVectorFromResults(Vec vector,int offset,int interp); @@ -210,6 +210,6 @@ ElementVector* CreatePVectorDiagnosticHutter(void); ElementMatrix* CreateJacobianDiagnosticMacayeal(void); - void GetSolutionFromInputsDiagnosticHoriz(Vec solution); - void GetSolutionFromInputsDiagnosticHutter(Vec solution); + void GetSolutionFromInputsDiagnosticHoriz(Vector* solution); + void GetSolutionFromInputsDiagnosticHutter(Vector* solution); void InputUpdateFromSolutionDiagnosticHoriz( double* solution); void InputUpdateFromSolutionDiagnosticHutter( double* solution); @@ -230,5 +230,5 @@ ElementVector* CreatePVectorHydrology(void); void CreateHydrologyWaterVelocityInput(void); - void GetSolutionFromInputsHydrology(Vec solution); + void GetSolutionFromInputsHydrology(Vector* solution); void InputUpdateFromSolutionHydrology(double* solution); #endif Index: /issm/trunk-jpl/src/c/objects/Loads/Icefront.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Loads/Icefront.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/objects/Loads/Icefront.cpp (revision 11679) @@ -316,5 +316,5 @@ /*}}}*/ /*FUNCTION Icefront::CreateKMatrix {{{1*/ -void Icefront::CreateKMatrix(Mat Kff, Mat Kfs){ +void Icefront::CreateKMatrix(Matrix* Kff, Matrix* Kfs){ /*No stiffness loads applied, do nothing: */ @@ -324,5 +324,5 @@ /*}}}*/ /*FUNCTION Icefront::CreatePVector {{{1*/ -void Icefront::CreatePVector(Vec pf){ +void Icefront::CreatePVector(Vector* pf){ /*Checks in debugging mode*/ @@ -362,10 +362,10 @@ /*}}}*/ /*FUNCTION Icefront::CreateJacobianMatrix{{{1*/ -void Icefront::CreateJacobianMatrix(Mat Jff){ +void Icefront::CreateJacobianMatrix(Matrix* Jff){ this->CreateKMatrix(Jff,NULL); } /*}}}1*/ /*FUNCTION Icefront::PenaltyCreateKMatrix {{{1*/ -void Icefront::PenaltyCreateKMatrix(Mat Kff, Mat Kfs, double kmax){ +void Icefront::PenaltyCreateKMatrix(Matrix* Kff, Matrix* Kfs, double kmax){ /*do nothing: */ return; @@ -373,5 +373,5 @@ /*}}}*/ /*FUNCTION Icefront::PenaltyCreatePVector{{{1*/ -void Icefront::PenaltyCreatePVector(Vec pf,double kmax){ +void Icefront::PenaltyCreatePVector(Vector* pf,double kmax){ /*do nothing: */ return; @@ -379,5 +379,5 @@ /*}}}*/ /*FUNCTION Icefront::PenaltyCreateJacobianMatrix{{{1*/ -void Icefront::PenaltyCreateJacobianMatrix(Mat Jff,double kmax){ +void Icefront::PenaltyCreateJacobianMatrix(Matrix* Jff,double kmax){ this->PenaltyCreateKMatrix(Jff,NULL,kmax); } Index: /issm/trunk-jpl/src/c/objects/Loads/Icefront.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Loads/Icefront.h (revision 11678) +++ /issm/trunk-jpl/src/c/objects/Loads/Icefront.h (revision 11679) @@ -74,10 +74,10 @@ void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); - void CreateKMatrix(Mat Kff, Mat Kfs); - void CreatePVector(Vec pf); - void CreateJacobianMatrix(Mat Jff); - void PenaltyCreateKMatrix(Mat Kff, Mat kfs, double kmax); - void PenaltyCreatePVector(Vec pf, double kmax); - void PenaltyCreateJacobianMatrix(Mat Jff,double kmax); + void CreateKMatrix(Matrix* Kff, Matrix* Kfs); + void CreatePVector(Vector* pf); + void CreateJacobianMatrix(Matrix* Jff); + void PenaltyCreateKMatrix(Matrix* Kff, Matrix* kfs, double kmax); + void PenaltyCreatePVector(Vector* pf, double kmax); + void PenaltyCreateJacobianMatrix(Matrix* Jff,double kmax); bool InAnalysis(int analysis_type); /*}}}*/ Index: /issm/trunk-jpl/src/c/objects/Loads/Load.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Loads/Load.h (revision 11678) +++ /issm/trunk-jpl/src/c/objects/Loads/Load.h (revision 11679) @@ -12,4 +12,6 @@ /*{{{1*/ class Object; +class Matrix; +class Vector; #include "../Object.h" @@ -27,10 +29,10 @@ virtual void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0; virtual void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0; - virtual void CreateKMatrix(Mat Kff, Mat Kfs)=0; - virtual void CreatePVector(Vec pf)=0; - virtual void CreateJacobianMatrix(Mat Jff)=0; - virtual void PenaltyCreateJacobianMatrix(Mat Jff,double kmax)=0; - virtual void PenaltyCreateKMatrix(Mat Kff, Mat Kfs, double kmax)=0; - virtual void PenaltyCreatePVector(Vec pf, double kmax)=0; + virtual void CreateKMatrix(Matrix* Kff, Matrix* Kfs)=0; + virtual void CreatePVector(Vector* pf)=0; + virtual void CreateJacobianMatrix(Matrix* Jff)=0; + virtual void PenaltyCreateJacobianMatrix(Matrix* Jff,double kmax)=0; + virtual void PenaltyCreateKMatrix(Matrix* Kff, Matrix* Kfs, double kmax)=0; + virtual void PenaltyCreatePVector(Vector* pf, double kmax)=0; virtual bool InAnalysis(int analysis_type)=0; /*}}}*/ Index: /issm/trunk-jpl/src/c/objects/Loads/Numericalflux.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Loads/Numericalflux.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/objects/Loads/Numericalflux.cpp (revision 11679) @@ -334,5 +334,5 @@ /*}}}*/ /*FUNCTION Numericalflux::CreateKMatrix {{{1*/ -void Numericalflux::CreateKMatrix(Mat Kff, Mat Kfs){ +void Numericalflux::CreateKMatrix(Matrix* Kff, Matrix* Kfs){ /*recover some parameters*/ @@ -365,5 +365,5 @@ /*}}}*/ /*FUNCTION Numericalflux::CreatePVector {{{1*/ -void Numericalflux::CreatePVector(Vec pf){ +void Numericalflux::CreatePVector(Vector* pf){ /*recover some parameters*/ @@ -395,5 +395,5 @@ /*}}}*/ /*FUNCTION Numericalflux::PenaltyCreateKMatrix {{{1*/ -void Numericalflux::PenaltyCreateKMatrix(Mat Kff, Mat Kfs,double kmax){ +void Numericalflux::PenaltyCreateKMatrix(Matrix* Kff, Matrix* Kfs,double kmax){ /*No stiffness loads applied, do nothing: */ @@ -403,5 +403,5 @@ /*}}}*/ /*FUNCTION Numericalflux::PenaltyCreatePVector{{{1*/ -void Numericalflux::PenaltyCreatePVector(Vec pf,double kmax){ +void Numericalflux::PenaltyCreatePVector(Vector* pf,double kmax){ /*No penalty loads applied, do nothing: */ Index: /issm/trunk-jpl/src/c/objects/Loads/Numericalflux.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Loads/Numericalflux.h (revision 11678) +++ /issm/trunk-jpl/src/c/objects/Loads/Numericalflux.h (revision 11679) @@ -70,10 +70,10 @@ void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); - void CreateKMatrix(Mat Kff, Mat Kfs); - void CreatePVector(Vec pf); - void CreateJacobianMatrix(Mat Jff){_error_("Not implemented yet");}; - void PenaltyCreateJacobianMatrix(Mat Jff,double kmax){_error_("Not implemented yet");}; - void PenaltyCreateKMatrix(Mat Kff, Mat kfs, double kmax); - void PenaltyCreatePVector(Vec pf, double kmax); + void CreateKMatrix(Matrix* Kff, Matrix* Kfs); + void CreatePVector(Vector* pf); + void CreateJacobianMatrix(Matrix* Jff){_error_("Not implemented yet");}; + void PenaltyCreateJacobianMatrix(Matrix* Jff,double kmax){_error_("Not implemented yet");}; + void PenaltyCreateKMatrix(Matrix* Kff, Matrix* kfs, double kmax); + void PenaltyCreatePVector(Vector* pf, double kmax); bool InAnalysis(int analysis_type); /*}}}*/ Index: /issm/trunk-jpl/src/c/objects/Loads/Pengrid.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Loads/Pengrid.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/objects/Loads/Pengrid.cpp (revision 11679) @@ -295,5 +295,5 @@ /*}}}1*/ /*FUNCTION Pengrid::CreateKMatrix {{{1*/ -void Pengrid::CreateKMatrix(Mat Kff, Mat Kfs){ +void Pengrid::CreateKMatrix(Matrix* Kff, Matrix* Kfs){ /*No loads applied, do nothing: */ @@ -303,5 +303,5 @@ /*}}}1*/ /*FUNCTION Pengrid::CreatePVector {{{1*/ -void Pengrid::CreatePVector(Vec pf){ +void Pengrid::CreatePVector(Vector* pf){ /*No loads applied, do nothing: */ @@ -311,5 +311,5 @@ /*}}}1*/ /*FUNCTION Pengrid::PenaltyCreateMatrix {{{1*/ -void Pengrid::PenaltyCreateKMatrix(Mat Kff, Mat Kfs,double kmax){ +void Pengrid::PenaltyCreateKMatrix(Matrix* Kff, Matrix* Kfs,double kmax){ /*Retrieve parameters: */ @@ -344,5 +344,5 @@ /*}}}1*/ /*FUNCTION Pengrid::PenaltyCreatePVector {{{1*/ -void Pengrid::PenaltyCreatePVector(Vec pf,double kmax){ +void Pengrid::PenaltyCreatePVector(Vector* pf,double kmax){ /*Retrieve parameters: */ Index: /issm/trunk-jpl/src/c/objects/Loads/Pengrid.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Loads/Pengrid.h (revision 11678) +++ /issm/trunk-jpl/src/c/objects/Loads/Pengrid.h (revision 11679) @@ -75,10 +75,10 @@ void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); - void CreateKMatrix(Mat Kff, Mat Kfs); - void CreatePVector(Vec pf); - void CreateJacobianMatrix(Mat Jff){_error_("Not implemented yet");}; - void PenaltyCreateJacobianMatrix(Mat Jff,double kmax){_error_("Not implemented yet");}; - void PenaltyCreateKMatrix(Mat Kff, Mat kfs, double kmax); - void PenaltyCreatePVector(Vec pf, double kmax); + void CreateKMatrix(Matrix* Kff, Matrix* Kfs); + void CreatePVector(Vector* pf); + void CreateJacobianMatrix(Matrix* Jff){_error_("Not implemented yet");}; + void PenaltyCreateJacobianMatrix(Matrix* Jff,double kmax){_error_("Not implemented yet");}; + void PenaltyCreateKMatrix(Matrix* Kff, Matrix* kfs, double kmax); + void PenaltyCreatePVector(Vector* pf, double kmax); bool InAnalysis(int analysis_type); /*}}}*/ Index: /issm/trunk-jpl/src/c/objects/Loads/Penpair.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Loads/Penpair.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/objects/Loads/Penpair.cpp (revision 11679) @@ -199,5 +199,5 @@ /*}}}1*/ /*FUNCTION Penpair::CreateKMatrix {{{1*/ -void Penpair::CreateKMatrix(Mat Kff, Mat Kfs){ +void Penpair::CreateKMatrix(Matrix* Kff, Matrix* Kfs){ /*If you code this piece, don't forget that a penalty will be inactive if it is dealing with clone nodes*/ /*No loads applied, do nothing: */ @@ -207,5 +207,5 @@ /*}}}1*/ /*FUNCTION Penpair::CreatePVector {{{1*/ -void Penpair::CreatePVector(Vec pf){ +void Penpair::CreatePVector(Vector* pf){ /*No loads applied, do nothing: */ @@ -215,10 +215,10 @@ /*}}}1*/ /*FUNCTION Penpair::CreateJacobianMatrix{{{1*/ -void Penpair::CreateJacobianMatrix(Mat Jff){ +void Penpair::CreateJacobianMatrix(Matrix* Jff){ this->CreateKMatrix(Jff,NULL); } /*}}}1*/ /*FUNCTION Penpair::PenaltyCreateKMatrix {{{1*/ -void Penpair::PenaltyCreateKMatrix(Mat Kff, Mat Kfs,double kmax){ +void Penpair::PenaltyCreateKMatrix(Matrix* Kff, Matrix* Kfs,double kmax){ /*Retrieve parameters: */ @@ -246,5 +246,5 @@ /*}}}1*/ /*FUNCTION Penpair::PenaltyCreatePVector {{{1*/ -void Penpair::PenaltyCreatePVector(Vec pf,double kmax){ +void Penpair::PenaltyCreatePVector(Vector* pf,double kmax){ /*No loads applied, do nothing: */ return; @@ -252,5 +252,5 @@ /*}}}1*/ /*FUNCTION Penpair::PenaltyCreateJacobianMatrix{{{1*/ -void Penpair::PenaltyCreateJacobianMatrix(Mat Jff,double kmax){ +void Penpair::PenaltyCreateJacobianMatrix(Matrix* Jff,double kmax){ this->PenaltyCreateKMatrix(Jff,NULL,kmax); } Index: /issm/trunk-jpl/src/c/objects/Loads/Penpair.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Loads/Penpair.h (revision 11678) +++ /issm/trunk-jpl/src/c/objects/Loads/Penpair.h (revision 11679) @@ -62,10 +62,10 @@ void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); - void CreateKMatrix(Mat Kff, Mat Kfs); - void CreatePVector(Vec pf); - void CreateJacobianMatrix(Mat Jff); - void PenaltyCreateKMatrix(Mat Kff,Mat Kfs,double kmax); - void PenaltyCreatePVector(Vec pf, double kmax); - void PenaltyCreateJacobianMatrix(Mat Jff,double kmax); + void CreateKMatrix(Matrix* Kff, Matrix* Kfs); + void CreatePVector(Vector* pf); + void CreateJacobianMatrix(Matrix* Jff); + void PenaltyCreateKMatrix(Matrix* Kff,Matrix* Kfs,double kmax); + void PenaltyCreatePVector(Vector* pf, double kmax); + void PenaltyCreateJacobianMatrix(Matrix* Jff,double kmax); bool InAnalysis(int analysis_type); /*}}}*/ Index: /issm/trunk-jpl/src/c/objects/Loads/Riftfront.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Loads/Riftfront.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/objects/Loads/Riftfront.cpp (revision 11679) @@ -428,5 +428,5 @@ /*}}}*/ /*FUNCTION Riftfront::PenaltyCreateKMatrix {{{1*/ -void Riftfront::PenaltyCreateKMatrix(Mat Kff, Mat Kfs,double kmax){ +void Riftfront::PenaltyCreateKMatrix(Matrix* Kff, Matrix* Kfs,double kmax){ /*Retrieve parameters: */ @@ -454,5 +454,5 @@ /*}}}1*/ /*FUNCTION Riftfront::PenaltyCreatePVector {{{1*/ -void Riftfront::PenaltyCreatePVector(Vec pf,double kmax){ +void Riftfront::PenaltyCreatePVector(Vector* pf,double kmax){ /*Retrieve parameters: */ @@ -480,5 +480,5 @@ /*}}}1*/ /*FUNCTION Riftfront::CreateKMatrix {{{1*/ -void Riftfront::CreateKMatrix(Mat Kff, Mat Kfs){ +void Riftfront::CreateKMatrix(Matrix* Kff, Matrix* Kfs){ /*do nothing: */ return; @@ -486,5 +486,5 @@ /*}}}1*/ /*FUNCTION Riftfront::CreatePVector {{{1*/ -void Riftfront::CreatePVector(Vec pf){ +void Riftfront::CreatePVector(Vector* pf){ /*do nothing: */ return; Index: /issm/trunk-jpl/src/c/objects/Loads/Riftfront.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Loads/Riftfront.h (revision 11678) +++ /issm/trunk-jpl/src/c/objects/Loads/Riftfront.h (revision 11679) @@ -82,10 +82,10 @@ void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); - void CreateKMatrix(Mat Kff, Mat Kfs); - void CreatePVector(Vec pf); - void CreateJacobianMatrix(Mat Jff){_error_("Not implemented yet");}; - void PenaltyCreateJacobianMatrix(Mat Jff,double kmax){_error_("Not implemented yet");}; - void PenaltyCreateKMatrix(Mat Kff, Mat kfs, double kmax); - void PenaltyCreatePVector(Vec pf, double kmax); + void CreateKMatrix(Matrix* Kff, Matrix* Kfs); + void CreatePVector(Vector* pf); + void CreateJacobianMatrix(Matrix* Jff){_error_("Not implemented yet");}; + void PenaltyCreateJacobianMatrix(Matrix* Jff,double kmax){_error_("Not implemented yet");}; + void PenaltyCreateKMatrix(Matrix* Kff, Matrix* kfs, double kmax); + void PenaltyCreatePVector(Vector* pf, double kmax); bool InAnalysis(int analysis_type); /*}}}*/ Index: /issm/trunk-jpl/src/c/objects/Node.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Node.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/objects/Node.cpp (revision 11679) @@ -594,5 +594,5 @@ /*}}}*/ /*FUNCTION Node::CreateNodalConstraints{{{1*/ -void Node::CreateNodalConstraints(Vec ys){ +void Node::CreateNodalConstraints(Vector* ys){ int i; @@ -613,5 +613,5 @@ /*Add values into constraint vector: */ - VecSetValues(ys,this->indexing.ssize,this->indexing.sdoflist,values,INSERT_VALUES); + ys->SetValues(this->indexing.ssize,this->indexing.sdoflist,values,INSERT_VALUES); } @@ -875,5 +875,5 @@ /*}}}*/ /*FUNCTION Node::VecMerge {{{1*/ -void Node::VecMerge(Vec ug, double* vector_serial,int setenum){ +void Node::VecMerge(Vector* ug, double* vector_serial,int setenum){ double* values=NULL; @@ -897,5 +897,5 @@ /*Add values into ug: */ - VecSetValues(ug,this->indexing.fsize,indices,(const double*)values,INSERT_VALUES); + ug->SetValues(this->indexing.fsize,indices,values,INSERT_VALUES); } } @@ -915,5 +915,5 @@ /*Add values into ug: */ - VecSetValues(ug,this->indexing.ssize,indices,(const double*)values,INSERT_VALUES); + ug->SetValues(this->indexing.ssize,indices,values,INSERT_VALUES); } } @@ -926,5 +926,5 @@ /*}}}*/ /*FUNCTION Node::VecReduce {{{1*/ -void Node::VecReduce(Vec vector, double* ug_serial,int setenum){ +void Node::VecReduce(Vector* vector, double* ug_serial,int setenum){ double* values=NULL; @@ -945,5 +945,5 @@ /*Add values into ug: */ - VecSetValues(vector,this->indexing.fsize,this->indexing.fdoflist,(const double*)values,INSERT_VALUES); + vector->SetValues(this->indexing.fsize,this->indexing.fdoflist,values,INSERT_VALUES); } } @@ -961,5 +961,5 @@ /*Add values into ug: */ - VecSetValues(vector,this->indexing.ssize,this->indexing.sdoflist,(const double*)values,INSERT_VALUES); + vector->SetValues(this->indexing.ssize,this->indexing.sdoflist,values,INSERT_VALUES); } } Index: /issm/trunk-jpl/src/c/objects/Node.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Node.h (revision 11678) +++ /issm/trunk-jpl/src/c/objects/Node.h (revision 11679) @@ -16,4 +16,6 @@ class DataSet; class Vertices; +class Vector; +class Matrix; #include "./Update.h" /*}}}*/ @@ -67,5 +69,5 @@ /*Node numerical routines {{{1*/ void Configure(DataSet* nodes,Vertices* vertices); - void CreateNodalConstraints(Vec ys); + void CreateNodalConstraints(Vector* ys); void SetCurrentConfiguration(DataSet* nodes,Vertices* vertices); int Sid(void); @@ -101,6 +103,6 @@ int IsGrounded(); void UpdateSpcs(double* ys); - void VecMerge(Vec ug, double* vector_serial,int setnum); - void VecReduce(Vec vector, double* ug_serial,int setnum); + void VecMerge(Vector* ug, double* vector_serial,int setenum); + void VecReduce(Vector* vector, double* ug_serial,int setnum); /*}}}*/ Index: /issm/trunk-jpl/src/c/objects/Numerics/ElementMatrix.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Numerics/ElementMatrix.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/objects/Numerics/ElementMatrix.cpp (revision 11679) @@ -245,6 +245,6 @@ /*ElementMatrix specific routines: */ -/*FUNCTION ElementMatrix::AddToGlobal(Mat Kff, Mat Kfs){{{1*/ -void ElementMatrix::AddToGlobal(Mat Kff, Mat Kfs){ +/*FUNCTION ElementMatrix::AddToGlobal(Matrix* Kff, Matrix* Kfs){{{1*/ +void ElementMatrix::AddToGlobal(Matrix* Kff, Matrix* Kfs){ int i,j; @@ -260,46 +260,50 @@ this->CheckConsistency(); - if(this->dofsymmetrical){ - /*only use row dofs to add values into global matrices: */ - - if(this->row_fsize){ - /*first, retrieve values that are in the f-set from the g-set values matrix: */ - localvalues=(double*)xmalloc(this->row_fsize*this->row_fsize*sizeof(double)); - for(i=0;irow_fsize;i++){ - for(j=0;jrow_fsize;j++){ - *(localvalues+this->row_fsize*i+j)=*(this->values+this->ncols*this->row_flocaldoflist[i]+this->row_flocaldoflist[j]); + #ifdef _HAVE_PETSC_ + if(this->dofsymmetrical){ + /*only use row dofs to add values into global matrices: */ + + if(this->row_fsize){ + /*first, retrieve values that are in the f-set from the g-set values matrix: */ + localvalues=(double*)xmalloc(this->row_fsize*this->row_fsize*sizeof(double)); + for(i=0;irow_fsize;i++){ + for(j=0;jrow_fsize;j++){ + *(localvalues+this->row_fsize*i+j)=*(this->values+this->ncols*this->row_flocaldoflist[i]+this->row_flocaldoflist[j]); + } } + /*add local values into global matrix, using the fglobaldoflist: */ + MatSetValues(Kff->matrix,this->row_fsize,this->row_fglobaldoflist,this->row_fsize,this->row_fglobaldoflist,(const double*)localvalues,ADD_VALUES); + + /*Free ressources:*/ + xfree((void**)&localvalues); } - /*add local values into global matrix, using the fglobaldoflist: */ - MatSetValues(Kff,this->row_fsize,this->row_fglobaldoflist,this->row_fsize,this->row_fglobaldoflist,(const double*)localvalues,ADD_VALUES); - - /*Free ressources:*/ - xfree((void**)&localvalues); - } - - - if((this->row_ssize!=0) && (this->row_fsize!=0)){ - /*first, retrieve values that are in the f and s-set from the g-set values matrix: */ - localvalues=(double*)xmalloc(this->row_fsize*this->row_ssize*sizeof(double)); - for(i=0;irow_fsize;i++){ - for(j=0;jrow_ssize;j++){ - *(localvalues+this->row_ssize*i+j)=*(this->values+this->ncols*this->row_flocaldoflist[i]+this->row_slocaldoflist[j]); + + + if((this->row_ssize!=0) && (this->row_fsize!=0)){ + /*first, retrieve values that are in the f and s-set from the g-set values matrix: */ + localvalues=(double*)xmalloc(this->row_fsize*this->row_ssize*sizeof(double)); + for(i=0;irow_fsize;i++){ + for(j=0;jrow_ssize;j++){ + *(localvalues+this->row_ssize*i+j)=*(this->values+this->ncols*this->row_flocaldoflist[i]+this->row_slocaldoflist[j]); + } } + /*add local values into global matrix, using the fglobaldoflist: */ + MatSetValues(Kfs->matrix,this->row_fsize,this->row_fglobaldoflist,this->row_ssize,this->row_sglobaldoflist,(const double*)localvalues,ADD_VALUES); + + /*Free ressources:*/ + xfree((void**)&localvalues); } - /*add local values into global matrix, using the fglobaldoflist: */ - MatSetValues(Kfs,this->row_fsize,this->row_fglobaldoflist,this->row_ssize,this->row_sglobaldoflist,(const double*)localvalues,ADD_VALUES); - - /*Free ressources:*/ - xfree((void**)&localvalues); - } - } - else{ - _error_(" non dofsymmetrical matrix AddToGlobal routine not support yet!"); - } - -} -/*}}}*/ -/*FUNCTION ElementMatrix::AddToGlobal(Mat Jff){{{1*/ -void ElementMatrix::AddToGlobal(Mat Jff){ + } + else{ + _error_(" non dofsymmetrical matrix AddToGlobal routine not support yet!"); + } + #else + _error_("not supported yet!"); + #endif + +} +/*}}}*/ +/*FUNCTION ElementMatrix::AddToGlobal(Matrix* Jff){{{1*/ +void ElementMatrix::AddToGlobal(Matrix* Jff){ int i,j; @@ -312,26 +316,30 @@ this->CheckConsistency(); - if(this->dofsymmetrical){ - /*only use row dofs to add values into global matrices: */ - - if(this->row_fsize){ - /*first, retrieve values that are in the f-set from the g-set values matrix: */ - localvalues=(double*)xmalloc(this->row_fsize*this->row_fsize*sizeof(double)); - for(i=0;irow_fsize;i++){ - for(j=0;jrow_fsize;j++){ - *(localvalues+this->row_fsize*i+j)=*(this->values+this->ncols*this->row_flocaldoflist[i]+this->row_flocaldoflist[j]); + #ifdef _HAVE_PETSC_ + if(this->dofsymmetrical){ + /*only use row dofs to add values into global matrices: */ + + if(this->row_fsize){ + /*first, retrieve values that are in the f-set from the g-set values matrix: */ + localvalues=(double*)xmalloc(this->row_fsize*this->row_fsize*sizeof(double)); + for(i=0;irow_fsize;i++){ + for(j=0;jrow_fsize;j++){ + *(localvalues+this->row_fsize*i+j)=*(this->values+this->ncols*this->row_flocaldoflist[i]+this->row_flocaldoflist[j]); + } } + /*add local values into global matrix, using the fglobaldoflist: */ + MatSetValues(Jff->matrix,this->row_fsize,this->row_fglobaldoflist,this->row_fsize,this->row_fglobaldoflist,(const double*)localvalues,ADD_VALUES); + + /*Free ressources:*/ + xfree((void**)&localvalues); } - /*add local values into global matrix, using the fglobaldoflist: */ - MatSetValues(Jff,this->row_fsize,this->row_fglobaldoflist,this->row_fsize,this->row_fglobaldoflist,(const double*)localvalues,ADD_VALUES); - - /*Free ressources:*/ - xfree((void**)&localvalues); - } - - } - else{ - _error_(" non dofsymmetrical matrix AddToGlobal routine not support yet!"); - } + + } + else{ + _error_(" non dofsymmetrical matrix AddToGlobal routine not support yet!"); + } + #else + _error_("not supported yet!"); + #endif } Index: /issm/trunk-jpl/src/c/objects/Numerics/ElementMatrix.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Numerics/ElementMatrix.h (revision 11678) +++ /issm/trunk-jpl/src/c/objects/Numerics/ElementMatrix.h (revision 11679) @@ -58,6 +58,6 @@ /*}}}*/ /*ElementMatrix specific routines {{{1*/ - void AddToGlobal(Mat Kff, Mat Kfs); - void AddToGlobal(Mat Jff); + void AddToGlobal(Matrix* Kff, Matrix* Kfs); + void AddToGlobal(Matrix* Jff); void Echo(void); void CheckConsistency(void); Index: /issm/trunk-jpl/src/c/objects/Numerics/ElementVector.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Numerics/ElementVector.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/objects/Numerics/ElementVector.cpp (revision 11679) @@ -160,42 +160,52 @@ /*ElementVector specific routines: */ -/*FUNCTION ElementVector::AddToGlobal(Vec pf){{{1*/ -void ElementVector::AddToGlobal(Vec pf){ +/*FUNCTION ElementVector::AddToGlobal(Vector* pf){{{1*/ +void ElementVector::AddToGlobal(Vector* pf){ int i; double* localvalues=NULL; - if(this->fsize){ - /*first, retrieve values that are in the f-set from the g-set values vector: */ - localvalues=(double*)xmalloc(this->fsize*sizeof(double)); - for(i=0;ifsize;i++){ - localvalues[i]=this->values[this->flocaldoflist[i]]; - } - /*add local values into global vector, using the fglobaldoflist: */ - VecSetValues(pf,this->fsize,this->fglobaldoflist,(const double*)localvalues,ADD_VALUES); - - /*Free ressources:*/ - xfree((void**)&localvalues); - } -} -/*}}}*/ -/*FUNCTION ElementVector::InsertIntoGlobal(Vec pf){{{1*/ -void ElementVector::InsertIntoGlobal(Vec pf){ + #ifdef _HAVE_PETSC_ + if(this->fsize){ + /*first, retrieve values that are in the f-set from the g-set values vector: */ + localvalues=(double*)xmalloc(this->fsize*sizeof(double)); + for(i=0;ifsize;i++){ + localvalues[i]=this->values[this->flocaldoflist[i]]; + } + /*add local values into global vector, using the fglobaldoflist: */ + VecSetValues(pf->vector,this->fsize,this->fglobaldoflist,(const double*)localvalues,ADD_VALUES); + + /*Free ressources:*/ + xfree((void**)&localvalues); + } + #else + _error_("not supported yet!"); + #endif + +} +/*}}}*/ +/*FUNCTION ElementVector::InsertIntoGlobal(Vector* pf){{{1*/ +void ElementVector::InsertIntoGlobal(Vector* pf){ int i; double* localvalues=NULL; - if(this->fsize){ - /*first, retrieve values that are in the f-set from the g-set values vector: */ - localvalues=(double*)xmalloc(this->fsize*sizeof(double)); - for(i=0;ifsize;i++){ - localvalues[i]=this->values[this->flocaldoflist[i]]; - } - /*add local values into global vector, using the fglobaldoflist: */ - VecSetValues(pf,this->fsize,this->fglobaldoflist,(const double*)localvalues,INSERT_VALUES); - - /*Free ressources:*/ - xfree((void**)&localvalues); - } + #ifdef _HAVE_PETSC_ + if(this->fsize){ + /*first, retrieve values that are in the f-set from the g-set values vector: */ + localvalues=(double*)xmalloc(this->fsize*sizeof(double)); + for(i=0;ifsize;i++){ + localvalues[i]=this->values[this->flocaldoflist[i]]; + } + /*add local values into global vector, using the fglobaldoflist: */ + VecSetValues(pf->vector,this->fsize,this->fglobaldoflist,(const double*)localvalues,INSERT_VALUES); + + /*Free ressources:*/ + xfree((void**)&localvalues); + } + #else + _error_("not supported yet!"); + #endif + } /*}}}*/ Index: /issm/trunk-jpl/src/c/objects/Numerics/ElementVector.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Numerics/ElementVector.h (revision 11678) +++ /issm/trunk-jpl/src/c/objects/Numerics/ElementVector.h (revision 11679) @@ -40,6 +40,6 @@ /*}}}*/ /*ElementVector specific routines {{{1*/ - void AddToGlobal(Vec pf); - void InsertIntoGlobal(Vec pf); + void AddToGlobal(Vector* pf); + void InsertIntoGlobal(Vector* pf); void Echo(void); void Init(ElementVector* pe); Index: /issm/trunk-jpl/src/c/objects/Numerics/Matrix.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Numerics/Matrix.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/objects/Numerics/Matrix.cpp (revision 11679) @@ -140,6 +140,63 @@ _error_("not implemented yet!"); #endif + return dataref; } /*}}}*/ #endif +/*FUNCTION Matrix::Assemble{{{1*/ +void Matrix::Assemble(void){ + #ifdef _HAVE_PETSC_ + MatAssemblyBegin(this->matrix,MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(this->matrix,MAT_FINAL_ASSEMBLY); + #if _PETSC_MAJOR_ == 2 + MatCompress(this->matrix); + #endif + #else + /*do nothing:*/ + #endif + +} +/*}}}*/ +/*FUNCTION Matrix::Norm{{{1*/ +double Matrix::Norm(int norm_type){ + + double norm=0; + #ifdef _HAVE_PETSC_ + MatNorm(this->matrix,(NormType)norm_type,&norm); + #else + _error_("not implemented yet!"); + #endif + return norm; +} +/*}}}*/ +/*FUNCTION Matrix::GetSize{{{1*/ +void Matrix::GetSize(int* pM,int* pN){ + + #ifdef _HAVE_PETSC_ + MatGetSize(this->matrix,pM,pN); + #else + _error_("not implemented yet!"); + #endif +} +/*}}}*/ +/*FUNCTION Matrix::GetLocalSize{{{1*/ +void Matrix::GetLocalSize(int* pM,int* pN){ + + #ifdef _HAVE_PETSC_ + MatGetLocalSize(this->matrix,pM,pN); + #else + _error_("not implemented yet!"); + #endif +} +/*}}}*/ +/*FUNCTION Matrix::MatMult{{{1*/ +void Matrix::MatMult(Vector* X,Vector* AX){ + + #ifdef _HAVE_PETSC_ + MatMultPatch(this->matrix,X->vector,AX->vector); + #else + _error_("not implemented yet!"); + #endif +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/objects/Numerics/Matrix.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Numerics/Matrix.h (revision 11678) +++ /issm/trunk-jpl/src/c/objects/Numerics/Matrix.h (revision 11679) @@ -25,4 +25,5 @@ #include "mex.h" #endif +class Vector; /*}}}*/ @@ -55,4 +56,9 @@ mxArray* ToMatlabMatrix(void); #endif + void Assemble(void); + double Norm(int norm_type); + void GetSize(int* pM,int* pN); + void GetLocalSize(int* pM,int* pN); + void MatMult(Vector* X,Vector* AX); /*}}}*/ Index: /issm/trunk-jpl/src/c/objects/Numerics/Vector.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Numerics/Vector.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/objects/Numerics/Vector.cpp (revision 11679) @@ -24,9 +24,9 @@ Vector::Vector(){ + this->M=0; #ifdef _HAVE_PETSC_ this->vector=NULL; #else this->vector=NULL; - this->M=0; #endif #ifdef _HAVE_ADOLC_ @@ -98,6 +98,151 @@ _error_("not implemented yet!"); #endif + return dataref; } /*}}}*/ #endif +/*FUNCTION Vector::Assemble{{{1*/ +void Vector::Assemble(void){ + + #ifdef _HAVE_PETSC_ + VecAssemblyBegin(this->vector); + VecAssemblyEnd(this->vector); + #else + /*do nothing*/ + #endif + +} +/*}}}*/ +/*FUNCTION Vector::SetValues{{{1*/ +void Vector::SetValues(int ssize, int* list, double* values, int mode){ + + + #ifdef _HAVE_PETSC_ + VecSetValues(this->vector,ssize,list,values,(InsertMode)mode); + #else + _error_("not implemented yet!"); + #endif + +} +/*}}}*/ +/*FUNCTION Vector::GetSize{{{1*/ +void Vector::GetSize(int* pM){ + + #ifdef _HAVE_PETSC_ + VecGetSize(this->vector,pM); + #else + _error_("not implemented yet!"); + #endif + +} +/*}}}*/ +/*FUNCTION Vector::GetLocalSize{{{1*/ +void Vector::GetLocalSize(int* pM){ + + #ifdef _HAVE_PETSC_ + VecGetLocalSize(this->vector,pM); + #else + _error_("not implemented yet!"); + #endif + +} +/*}}}*/ +/*FUNCTION Vector::Duplicate{{{1*/ +Vector* Vector::Duplicate(void){ + + Vector* output=NULL; + output=new Vector(); + #ifdef _HAVE_PETSC_ + Vec vec_output=NULL; + VecDuplicate(this->vector,&vec_output); + output->vector=vec_output; + VecGetSize(output->vector,&output->M); + #else + _error_("not implemented yet!"); + #endif + return output; + +} +/*}}}*/ +/*FUNCTION Vector::Set{{{1*/ +void Vector::Set(double value){ + + #ifdef _HAVE_PETSC_ + this->Set(value); + #else + _error_("not implemented yet!"); + #endif + +} +/*}}}*/ +/*FUNCTION Vector::AXPY{{{1*/ +void Vector::AXPY(Vector* X, double a){ + + #ifdef _HAVE_PETSC_ + VecAXPY(this->vector,a,X->vector); + #else + _error_("not implemented yet!"); + #endif + +} +/*}}}*/ +/*FUNCTION Vector::AYPX{{{1*/ +void Vector::AYPX(Vector* X, double a){ + + #ifdef _HAVE_PETSC_ + VecAYPX(this->vector,a,X->vector); + #else + _error_("not implemented yet!"); + #endif + +} +/*}}}*/ +/*FUNCTION Vector::ToMPISerial{{{1*/ +double* Vector::ToMPISerial(void){ + + double* vec_serial=NULL; + + #ifdef _HAVE_PETSC_ + VecToMPISerial(&vec_serial, this->vector); + #else + _error_("not implemented yet!"); + #endif + + return vec_serial; + +} +/*}}}*/ +/*FUNCTION Vector::Copy{{{1*/ +void Vector::Copy(Vector* to){ + + #ifdef _HAVE_PETSC_ + VecCopy(this->vector,to->vector); + #else + _error_("not implemented yet!"); + #endif + +} +/*}}}*/ +/*FUNCTION Vector::Norm{{{1*/ +double Vector::Norm(int norm_type){ + + double norm=0; + #ifdef _HAVE_PETSC_ + VecNorm(this->vector,(NormType)norm_type,&norm); + #else + _error_("not implemented yet!"); + #endif + return norm; +} +/*}}}*/ +/*FUNCTION Vector::Scale{{{1*/ +void Vector::Scale(double scale_factor){ + + #ifdef _HAVE_PETSC_ + VecScale(this->vector,scale_factor); + #else + _error_("not implemented yet!"); + #endif +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/objects/Numerics/Vector.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Numerics/Vector.h (revision 11678) +++ /issm/trunk-jpl/src/c/objects/Numerics/Vector.h (revision 11679) @@ -54,4 +54,16 @@ mxArray* ToMatlabVector(void); #endif + void Assemble(void); + void SetValues(int ssize, int* list, double* values, int mode); + void GetSize(int* pM); + void GetLocalSize(int* pM); + Vector* Duplicate(void); + void Set(double value); + void AXPY(Vector* X, double a); + void AYPX(Vector* X, double a); + double* ToMPISerial(void); + void Copy(Vector* to); + double Norm(int norm_type); + void Scale(double scale_factor); /*}}}*/ }; Index: /issm/trunk-jpl/src/c/shared/Alloc/alloc.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Alloc/alloc.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/shared/Alloc/alloc.cpp (revision 11679) @@ -29,4 +29,5 @@ #include "../Exceptions/exceptions.h" #include "../../include/include.h" +#include "../../objects/objects.h" void* xmalloc(int size){ @@ -80,4 +81,23 @@ } +void xdelete( Matrix** pv){ + + if (pv && *pv) { + + delete *pv; + *pv=NULL; + } +} + +void xdelete( Vector** pv){ + + if (pv && *pv) { + + delete *pv; + *pv=NULL; + } +} + + void* xrealloc ( void* pv, int size){ Index: /issm/trunk-jpl/src/c/shared/Alloc/alloc.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Alloc/alloc.h (revision 11678) +++ /issm/trunk-jpl/src/c/shared/Alloc/alloc.h (revision 11679) @@ -5,9 +5,12 @@ #ifndef _ALLOC_H_ #define _ALLOC_H_ - +class Matrix; +class Vector; void* xmalloc(int size); void* xcalloc(int n,int size); void xfree(void** pvptr); void* xrealloc ( void* pv, int size); +void xdelete(Matrix** pvptr); +void xdelete(Vector** pvptr); #endif Index: /issm/trunk-jpl/src/c/solutions/ResetBoundaryConditions.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/ResetBoundaryConditions.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/solutions/ResetBoundaryConditions.cpp (revision 11679) @@ -11,5 +11,5 @@ /*variables: */ - Vec yg = NULL; + Vector* yg = NULL; Nodes *nodes = NULL; int i; @@ -30,4 +30,4 @@ /*Free ressources:*/ - VecFree(&yg); + delete yg; } Index: /issm/trunk-jpl/src/c/solutions/convergence.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/convergence.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/solutions/convergence.cpp (revision 11679) @@ -8,5 +8,5 @@ #include "../EnumDefinitions/EnumDefinitions.h" -void convergence(bool* pconverged, Mat Kff,Vec pf,Vec uf,Vec old_uf,Parameters* parameters){ +void convergence(bool* pconverged, Matrix* Kff,Vector* pf,Vector* uf,Vector* old_uf,Parameters* parameters){ /*output*/ @@ -14,9 +14,9 @@ /*intermediary*/ - Vec KU=NULL; - Vec KUF=NULL; - Vec KUold=NULL; - Vec KUoldF=NULL; - Vec duf=NULL; + Vector* KU=NULL; + Vector* KUF=NULL; + Vector* KUold=NULL; + Vector* KUoldF=NULL; + Vector* duf=NULL; double ndu,nduinf,nu; double nKUF; @@ -48,16 +48,16 @@ //compute KUF = KU - F = K*U - F - VecDuplicate(uf,&KU); MatMultPatch(Kff,uf,KU); - VecDuplicate(KU,&KUF);VecCopy(KU,KUF); VecAYPX(KUF,-1.0,pf); + KU=uf->Duplicate(); Kff->MatMult(uf,KU); + KUF=KU->Duplicate(); KU->Copy(KUF); KUF->AYPX(pf,-1.0); //compute norm(KUF), norm(F) and residue - VecNorm(KUF,NORM_2,&nKUF); - VecNorm(pf,NORM_2,&nF); + nKUF=KUF->Norm(NORM_2); + nF=pf->Norm(NORM_2); solver_residue=nKUF/nF; _printf_(true,"\n%s%g\n"," solver residue: norm(KU-F)/norm(F)=",solver_residue); //clean up - VecFree(&KU); - VecFree(&KUF); + delete KU; + delete KUF; } @@ -65,8 +65,8 @@ //compute K[n]U[n-1]F = K[n]U[n-1] - F - VecDuplicate(uf,&KUold); MatMultPatch(Kff,old_uf,KUold); - VecDuplicate(KUold,&KUoldF);VecCopy(KUold,KUoldF); VecAYPX(KUoldF,-1.0,pf); - VecNorm(KUoldF,NORM_2,&nKUoldF); - VecNorm(pf,NORM_2,&nF); + KUold=uf->Duplicate(); Kff->MatMult(old_uf,KUold); + KUoldF=KUold->Duplicate();KUold->Copy(KUoldF); KUoldF->AYPX(pf,-1.0); + nKUoldF=KUoldF->Norm(NORM_2); + nF=pf->Norm(NORM_2); res=nKUoldF/nF; if (isnan(res)){ @@ -76,6 +76,6 @@ //clean up - VecFree(&KUold); - VecFree(&KUoldF); + delete KUold; + delete KUoldF; //print @@ -93,11 +93,11 @@ //compute norm(du)/norm(u) - VecDuplicate(old_uf,&duf);VecCopy(old_uf,duf); VecAYPX(duf,-1.0,uf); - VecNorm(duf,NORM_2,&ndu); VecNorm(old_uf,NORM_2,&nu); + duf=old_uf->Duplicate(); old_uf->Copy(duf); duf->AYPX(uf,-1.0); + ndu=duf->Norm(NORM_2); nu=old_uf->Norm(NORM_2); if (isnan(ndu) || isnan(nu)) _error_("convergence criterion is NaN!"); //clean up - VecFree(&duf); + delete duf; //print @@ -119,10 +119,10 @@ //compute max(du) - VecDuplicate(old_uf,&duf);VecCopy(old_uf,duf); VecAYPX(duf,-1.0,uf); - VecNorm(duf,NORM_2,&ndu); VecNorm(duf,NORM_INFINITY,&nduinf); + duf=old_uf->Duplicate(); old_uf->Copy(duf); duf->AYPX(uf,-1.0); + ndu=duf->Norm(NORM_2); nduinf=duf->Norm(NORM_INFINITY); if (isnan(ndu) || isnan(nu)) _error_("convergence criterion is NaN!"); //clean up - VecFree(&duf); + delete duf; //print Index: /issm/trunk-jpl/src/c/solutions/solutions.h =================================================================== --- /issm/trunk-jpl/src/c/solutions/solutions.h (revision 11678) +++ /issm/trunk-jpl/src/c/solutions/solutions.h (revision 11679) @@ -35,5 +35,5 @@ //convergence: -void convergence(bool* pconverged, Mat K_ff,Vec p_f,Vec u_f,Vec u_f_old,Parameters* parameters); +void convergence(bool* pconverged, Matrix* K_ff,Vector* p_f,Vector* u_f,Vector* u_f_old,Parameters* parameters); bool controlconvergence(double J,double tol_cm); bool steadystateconvergence(FemModel* femmodel); Index: /issm/trunk-jpl/src/c/solvers/solver_adjoint_linear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solvers/solver_adjoint_linear.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/solvers/solver_adjoint_linear.cpp (revision 11679) @@ -13,9 +13,11 @@ /*intermediary: */ - Mat Kff = NULL, Kfs = NULL; - Vec ug = NULL, uf = NULL; - Vec pf = NULL; - Vec df = NULL; - Vec ys = NULL; + Matrix* Kff = NULL; + Matrix* Kfs = NULL; + Vector* ug = NULL; + Vector* uf = NULL; + Vector* pf = NULL; + Vector* df = NULL; + Vector* ys = NULL; int configuration_type; @@ -26,8 +28,8 @@ 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,true); MatFree(&Kfs); //true means spc = 0 - Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); MatFree(&Kff); VecFree(&pf); VecFree(&df); - Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters,true); VecFree(&uf);VecFree(&ys); //true means spc0 + Reduceloadx(pf, Kfs, ys,true); xdelete(&Kfs); //true means spc = 0 + Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); xdelete(&Kff); xdelete(&pf); xdelete(&df); + Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters,true); xdelete(&uf);xdelete(&ys); //true means spc0 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); - VecFree(&ug); VecFree(&uf); + xdelete(&ug); xdelete(&uf); } Index: /issm/trunk-jpl/src/c/solvers/solver_linear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solvers/solver_linear.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/solvers/solver_linear.cpp (revision 11679) @@ -11,10 +11,11 @@ /*intermediary: */ - Mat Kff = NULL, Kfs = NULL; - Vec ug = NULL; - Vec uf = NULL; - Vec pf = NULL; - Vec df = NULL; - Vec ys = NULL; + Matrix* Kff = NULL; + Matrix* Kfs = NULL; + Vector* ug = NULL; + Vector* uf = NULL; + Vector* pf = NULL; + Vector* df = NULL; + Vector* ys = NULL; int configuration_type; @@ -25,8 +26,9 @@ 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, NULL, df, femmodel->parameters); MatFree(&Kff); VecFree(&pf); VecFree(&df); - Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);VecFree(&uf);VecFree(&ys); + Reduceloadx(pf, Kfs, ys); xdelete(&Kfs); + Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); + xdelete(&Kff); xdelete(&pf); xdelete(&df); + Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);xdelete(&uf); xdelete(&ys); InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); - VecFree(&ug); VecFree(&uf); + xdelete(&ug); xdelete(&uf); } Index: /issm/trunk-jpl/src/c/solvers/solver_newton.cpp =================================================================== --- /issm/trunk-jpl/src/c/solvers/solver_newton.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/solvers/solver_newton.cpp (revision 11679) @@ -18,10 +18,16 @@ int count; double kmax; - 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; + Matrix* Kff = NULL; + Matrix* Kfs = NULL; + Matrix* Jff = NULL; + Vector* ug = NULL; + Vector* old_ug = NULL; + Vector* uf = NULL; + Vector* old_uf = NULL; + Vector* duf = NULL; + Vector* pf = NULL; + Vector* pJf = NULL; + Vector* df = NULL; + Vector* ys = NULL; /*parameters:*/ @@ -47,18 +53,18 @@ for(;;){ - VecFree(&old_ug);old_ug=ug; - VecFree(&old_uf);old_uf=uf; + xdelete(&old_ug);old_ug=ug; + xdelete(&old_uf);old_uf=uf; /*Solver forward model*/ 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);VecFree(&df); - Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);VecFree(&ys); - InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);VecFree(&ug); + Reduceloadx(pf,Kfs,ys);xdelete(&Kfs); + Solverx(&uf,Kff,pf,old_uf,df,femmodel->parameters);xdelete(&df); + Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);xdelete(&ys); + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);xdelete(&ug); /*Check convergence*/ convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters); - MatFree(&Kff);VecFree(&pf); + xdelete(&Kff); xdelete(&pf); if(converged==true) break; if(count>=max_nonlinear_iterations){ @@ -70,15 +76,13 @@ SystemMatricesx(&Kff,&Kfs,&pf,NULL,&kmax,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); - Reduceloadx(pf,Kfs,ys); MatFree(&Kfs); + Reduceloadx(pf,Kfs,ys); xdelete(&Kfs); - VecDuplicate(pf,&pJf); - MatMultPatch(Kff,uf,pJf); MatFree(&Kff); - VecScale(pJf,-1.); - VecAXPY(pJf,+1.,pf); VecFree(&pf); + pJf=pf->Duplicate(); Kff->MatMult(uf,pJf); xdelete(&Kff); + pJf->Scale(-1.0); pJf->AXPY(pf,+1.0); xdelete(&pf); CreateJacobianMatrixx(&Jff,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,kmax); - Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters); MatFree(&Jff);VecFree(&pJf); - VecAXPY(uf,1.,duf); VecFree(&duf); - Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);VecFree(&ys); + Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters); xdelete(&Jff); xdelete(&pJf); + uf->AXPY(duf, 1.0); xdelete(&duf); + Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);xdelete(&ys); InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug); @@ -89,7 +93,7 @@ /*clean-up*/ - VecFree(&uf); - VecFree(&ug); - VecFree(&old_ug); - VecFree(&old_uf); + xdelete(&uf); + xdelete(&ug); + xdelete(&old_ug); + xdelete(&old_uf); } Index: /issm/trunk-jpl/src/c/solvers/solver_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solvers/solver_nonlinear.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/solvers/solver_nonlinear.cpp (revision 11679) @@ -14,10 +14,13 @@ /*intermediary: */ - Mat Kff = NULL, Kfs = NULL; - Vec ug = NULL, old_ug = NULL; - Vec uf = NULL, old_uf = NULL; - Vec pf = NULL; - Vec df = NULL; - Vec ys = NULL; + Matrix* Kff = NULL; + Matrix* Kfs = NULL; + Vector* ug = NULL; + Vector* old_ug = NULL; + Vector* uf = NULL; + Vector* old_uf = NULL; + Vector* pf = NULL; + Vector* df = NULL; + Vector* ys = NULL; Loads* loads=NULL; @@ -56,14 +59,14 @@ //save pointer to old velocity - VecFree(&old_ug);old_ug=ug; - VecFree(&old_uf);old_uf=uf; + xdelete(&old_ug);old_ug=ug; + xdelete(&old_uf);old_uf=uf; SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters); CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); - Reduceloadx(pf, Kfs, ys); MatFree(&Kfs); + Reduceloadx(pf, Kfs, ys); xdelete(&Kfs); Solverx(&uf, Kff, pf, old_uf, df, femmodel->parameters); - Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);VecFree(&ys); + Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);xdelete(&ys); - convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters); MatFree(&Kff);VecFree(&pf); VecFree(&df); + convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters); xdelete(&Kff); xdelete(&pf); xdelete(&df); 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); @@ -96,7 +99,7 @@ /*clean-up*/ if(conserve_loads) delete loads; - VecFree(&uf); - VecFree(&ug); - VecFree(&old_ug); - VecFree(&old_uf); + xdelete(&uf); + xdelete(&ug); + xdelete(&old_ug); + xdelete(&old_uf); } Index: /issm/trunk-jpl/src/c/solvers/solver_stokescoupling_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solvers/solver_stokescoupling_nonlinear.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/solvers/solver_stokescoupling_nonlinear.cpp (revision 11679) @@ -14,13 +14,18 @@ /*intermediary: */ - Mat Kff_horiz = NULL, Kfs_horiz = NULL; - Vec ug_horiz = NULL, uf_horiz = NULL, old_uf_horiz = NULL; - Vec pf_horiz = NULL; - Vec df_horiz = NULL; - Mat Kff_vert = NULL, Kfs_vert = NULL; - Vec ug_vert = NULL, uf_vert = NULL; - Vec pf_vert = NULL; - Vec df_vert = NULL; - Vec ys = NULL; + Matrix* Kff_horiz = NULL; + Matrix* Kfs_horiz = NULL; + Vector* ug_horiz = NULL; + Vector* uf_horiz = NULL; + Vector* old_uf_horiz = NULL; + Vector* pf_horiz = NULL; + Vector* df_horiz = NULL; + Matrix* Kff_vert = NULL; + Matrix* Kfs_vert = NULL; + Vector* ug_vert = NULL; + Vector* uf_vert = NULL; + Vector* pf_vert = NULL; + Vector* df_vert = NULL; + Vector* ys = NULL; bool converged; int constraints_converged; @@ -54,18 +59,18 @@ //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_horiz); - VecFree(&ug_horiz); + xdelete(&ug_horiz); //save pointer to old velocity - VecFree(&old_uf_horiz);old_uf_horiz=uf_horiz; + xdelete(&old_uf_horiz); old_uf_horiz=uf_horiz; /*solve: */ SystemMatricesx(&Kff_horiz, &Kfs_horiz, &pf_horiz, &df_horiz, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); - Reduceloadx(pf_horiz, Kfs_horiz, ys); MatFree(&Kfs_horiz); + Reduceloadx(pf_horiz, Kfs_horiz, ys); xdelete(&Kfs_horiz); Solverx(&uf_horiz, Kff_horiz, pf_horiz, old_uf_horiz, df_horiz,femmodel->parameters); - Mergesolutionfromftogx(&ug_horiz, uf_horiz,ys,femmodel->nodes,femmodel->parameters); VecFree(&ys); + Mergesolutionfromftogx(&ug_horiz, uf_horiz,ys,femmodel->nodes,femmodel->parameters); xdelete(&ys); InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug_horiz); - convergence(&converged,Kff_horiz,pf_horiz,uf_horiz,old_uf_horiz,femmodel->parameters); MatFree(&Kff_horiz);VecFree(&pf_horiz); VecFree(&df_horiz); + convergence(&converged,Kff_horiz,pf_horiz,uf_horiz,old_uf_horiz,femmodel->parameters); xdelete(&Kff_horiz); xdelete(&pf_horiz); xdelete(&df_horiz); /*Second compute vertical velocity: */ @@ -76,9 +81,9 @@ SystemMatricesx(&Kff_vert, &Kfs_vert, &pf_vert, &df_vert,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); - Reduceloadx(pf_vert, Kfs_vert, ys); MatFree(&Kfs_vert); - Solverx(&uf_vert, Kff_vert, pf_vert, NULL, df_vert,femmodel->parameters); MatFree(&Kff_vert); VecFree(&pf_vert); VecFree(&df_vert); - Mergesolutionfromftogx(&ug_vert, uf_vert,ys,femmodel->nodes,femmodel->parameters);VecFree(&uf_vert); VecFree(&ys); + Reduceloadx(pf_vert, Kfs_vert, ys); xdelete(&Kfs_vert); + Solverx(&uf_vert, Kff_vert, pf_vert, NULL, df_vert,femmodel->parameters); xdelete(&Kff_vert); xdelete(&pf_vert); xdelete(&df_vert); + Mergesolutionfromftogx(&ug_vert, uf_vert,ys,femmodel->nodes,femmodel->parameters);xdelete(&uf_vert); xdelete(&ys); InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug_vert); - VecFree(&ug_vert); VecFree(&uf_vert); + xdelete(&ug_vert); xdelete(&uf_vert); /*Increase count: */ @@ -92,7 +97,7 @@ /*clean-up*/ - VecFree(&old_uf_horiz); - VecFree(&uf_horiz); - VecFree(&ug_horiz); - VecFree(&ys); + xdelete(&old_uf_horiz); + xdelete(&uf_horiz); + xdelete(&ug_horiz); + xdelete(&ys); } Index: /issm/trunk-jpl/src/c/solvers/solver_thermal_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solvers/solver_thermal_nonlinear.cpp (revision 11678) +++ /issm/trunk-jpl/src/c/solvers/solver_thermal_nonlinear.cpp (revision 11679) @@ -12,15 +12,15 @@ /*solution : */ - Vec tg=NULL; - Vec tf=NULL; - Vec tf_old=NULL; - Vec ys=NULL; + Vector* tg=NULL; + Vector* tf=NULL; + Vector* tf_old=NULL; + Vector* ys=NULL; double melting_offset; /*intermediary: */ - Mat Kff=NULL; - Mat Kfs=NULL; - Vec pf=NULL; - Vec df=NULL; + Matrix* Kff=NULL; + Matrix* Kfs=NULL; + Vector* pf=NULL; + Vector* df=NULL; bool converged; @@ -56,9 +56,9 @@ SystemMatricesx(&Kff, &Kfs, &pf,&df, &melting_offset,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); - Reduceloadx(pf, Kfs, ys); MatFree(&Kfs); VecFree(&tf); + Reduceloadx(pf, Kfs, ys); xdelete(&Kfs); xdelete(&tf); Solverx(&tf, Kff, pf,tf_old, df, femmodel->parameters); - VecFree(&tf_old); VecDuplicatePatch(&tf_old,tf); - MatFree(&Kff);VecFree(&pf);VecFree(&tg); VecFree(&df); - Mergesolutionfromftogx(&tg, tf,ys,femmodel->nodes,femmodel->parameters); VecFree(&ys); + xdelete(&tf_old); tf_old=tf->Duplicate(); + xdelete(&Kff);xdelete(&pf);xdelete(&tg); xdelete(&df); + Mergesolutionfromftogx(&tg, tf,ys,femmodel->nodes,femmodel->parameters); xdelete(&ys); InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,tg); @@ -84,7 +84,7 @@ /*Free ressources: */ - VecFree(&tg); - VecFree(&tf); - VecFree(&tf_old); - VecFree(&ys); + xdelete(&tg); + xdelete(&tf); + xdelete(&tf_old); + xdelete(&ys); } Index: /issm/trunk-jpl/src/mex/CreateJacobianMatrix/CreateJacobianMatrix.cpp =================================================================== --- /issm/trunk-jpl/src/mex/CreateJacobianMatrix/CreateJacobianMatrix.cpp (revision 11678) +++ /issm/trunk-jpl/src/mex/CreateJacobianMatrix/CreateJacobianMatrix.cpp (revision 11679) @@ -17,5 +17,5 @@ /* output datasets: */ - Mat Jff = NULL; + Matrix* Jff = NULL; /*Boot module: */ @@ -53,5 +53,5 @@ delete materials; delete parameters; - MatFree(&Jff); + delete Jff; /*end module: */ Index: /issm/trunk-jpl/src/mex/CreateNodalConstraints/CreateNodalConstraints.cpp =================================================================== --- /issm/trunk-jpl/src/mex/CreateNodalConstraints/CreateNodalConstraints.cpp (revision 11678) +++ /issm/trunk-jpl/src/mex/CreateNodalConstraints/CreateNodalConstraints.cpp (revision 11679) @@ -12,5 +12,5 @@ /* output datasets: */ - Vec ys=NULL; + Vector* ys=NULL; /*Boot module: */ @@ -32,5 +32,5 @@ /*Free ressources: */ delete nodes; - VecFree(&ys); + delete ys; /*end module: */ Index: /issm/trunk-jpl/src/mex/GetSolutionFromInputs/GetSolutionFromInputs.cpp =================================================================== --- /issm/trunk-jpl/src/mex/GetSolutionFromInputs/GetSolutionFromInputs.cpp (revision 11678) +++ /issm/trunk-jpl/src/mex/GetSolutionFromInputs/GetSolutionFromInputs.cpp (revision 11679) @@ -14,5 +14,5 @@ Materials* materials=NULL; Parameters* parameters=NULL; - Vec ug=NULL; + Vector* ug=NULL; /* output datasets: elements and loads*/ Index: /issm/trunk-jpl/src/mex/InputUpdateFromSolution/InputUpdateFromSolution.cpp =================================================================== --- /issm/trunk-jpl/src/mex/InputUpdateFromSolution/InputUpdateFromSolution.cpp (revision 11678) +++ /issm/trunk-jpl/src/mex/InputUpdateFromSolution/InputUpdateFromSolution.cpp (revision 11679) @@ -14,5 +14,5 @@ Materials* materials=NULL; Parameters* parameters=NULL; - Vec solution=NULL; + Vector* solution=NULL; /*Boot module: */ @@ -50,5 +50,5 @@ delete materials; delete parameters; - VecFree(&solution); + delete solution; /*end module: */ Index: /issm/trunk-jpl/src/mex/Mergesolutionfromftog/Mergesolutionfromftog.cpp =================================================================== --- /issm/trunk-jpl/src/mex/Mergesolutionfromftog/Mergesolutionfromftog.cpp (revision 11678) +++ /issm/trunk-jpl/src/mex/Mergesolutionfromftog/Mergesolutionfromftog.cpp (revision 11679) @@ -9,11 +9,11 @@ /*input datasets: */ bool flag_ys0; - Vec uf = NULL; - Vec ys = NULL; + Vector* uf = NULL; + Vector* ys = NULL; Nodes* nodes = NULL; Parameters* parameters = NULL; /* output datasets: */ - Vec ug=NULL; + Vector* ug=NULL; /*Boot module: */ @@ -45,7 +45,7 @@ /*Free ressources: */ - VecFree(&uf); - VecFree(&ug); - VecFree(&ys); + delete uf; + delete ug; + delete ys; delete nodes; delete parameters; Index: /issm/trunk-jpl/src/mex/Reduceload/Reduceload.cpp =================================================================== --- /issm/trunk-jpl/src/mex/Reduceload/Reduceload.cpp (revision 11678) +++ /issm/trunk-jpl/src/mex/Reduceload/Reduceload.cpp (revision 11679) @@ -8,7 +8,7 @@ /*input datasets: */ - Vec pf = NULL; - Mat Kfs = NULL; - Vec ys = NULL; + Vector* pf = NULL; + Matrix* Kfs = NULL; + Vector* ys = NULL; bool flag_ys0=false; @@ -40,7 +40,7 @@ /*Free ressources: */ - VecFree(&pf); - MatFree(&Kfs); - VecFree(&ys); + delete pf; + delete Kfs; + delete ys; MODULEEND(); Index: /issm/trunk-jpl/src/mex/Reducevectorgtof/Reducevectorgtof.cpp =================================================================== --- /issm/trunk-jpl/src/mex/Reducevectorgtof/Reducevectorgtof.cpp (revision 11678) +++ /issm/trunk-jpl/src/mex/Reducevectorgtof/Reducevectorgtof.cpp (revision 11679) @@ -8,10 +8,10 @@ /*input datasets: */ - Vec ug=NULL; + Vector* ug=NULL; Nodes* nodes=NULL; Parameters* parameters=NULL; /* output datasets: */ - Vec uf=NULL; + Vector* uf=NULL; /*Boot module: */ @@ -35,6 +35,6 @@ delete nodes; delete parameters; - VecFree(&ug); - VecFree(&uf); + delete ug; + delete uf; /*end module: */ Index: /issm/trunk-jpl/src/mex/Reducevectorgtos/Reducevectorgtos.cpp =================================================================== --- /issm/trunk-jpl/src/mex/Reducevectorgtos/Reducevectorgtos.cpp (revision 11678) +++ /issm/trunk-jpl/src/mex/Reducevectorgtos/Reducevectorgtos.cpp (revision 11679) @@ -8,10 +8,10 @@ /*input datasets: */ - Vec yg=NULL; + Vector* yg=NULL; Nodes* nodes=NULL; Parameters* parameters=NULL; /* output datasets: */ - Vec ys=NULL; + Vector* ys=NULL; /*Boot module: */ @@ -35,6 +35,6 @@ delete nodes; delete parameters; - VecFree(&yg); - VecFree(&ys); + delete yg; + delete ys; /*end module: */ Index: /issm/trunk-jpl/src/mex/Solver/Solver.cpp =================================================================== --- /issm/trunk-jpl/src/mex/Solver/Solver.cpp (revision 11678) +++ /issm/trunk-jpl/src/mex/Solver/Solver.cpp (revision 11679) @@ -8,9 +8,9 @@ /*input datasets: */ - Mat Kff = NULL; - Vec pf = NULL; - Vec uf0 = NULL; - Vec uf = NULL; - Vec df = NULL; + Matrix* Kff = NULL; + Vector* pf = NULL; + Vector* uf0 = NULL; + Vector* uf = NULL; + Vector* df = NULL; Parameters *parameters = NULL; int analysis_type; @@ -66,9 +66,9 @@ /*Free ressources: */ - MatFree(&Kff); - VecFree(&pf); - VecFree(&uf0); - VecFree(&uf); - VecFree(&df); + delete Kff; + delete pf; + delete uf0; + delete uf; + delete df; delete parameters; Index: /issm/trunk-jpl/src/mex/SystemMatrices/SystemMatrices.cpp =================================================================== --- /issm/trunk-jpl/src/mex/SystemMatrices/SystemMatrices.cpp (revision 11678) +++ /issm/trunk-jpl/src/mex/SystemMatrices/SystemMatrices.cpp (revision 11679) @@ -17,8 +17,8 @@ /* output datasets: */ - Mat Kff = NULL; - Mat Kfs = NULL; - Vec pf = NULL; - Vec df = NULL; + Matrix* Kff = NULL; + Matrix* Kfs = NULL; + Vector* pf = NULL; + Vector* df = NULL; double kmax; @@ -72,8 +72,8 @@ delete materials; delete parameters; - MatFree(&Kff); - MatFree(&Kfs); - VecFree(&pf); - VecFree(&df); + delete Kff; + delete Kfs; + delete pf; + delete df; /*end module: */ Index: /issm/trunk-jpl/src/mex/UpdateDynamicConstraints/UpdateDynamicConstraints.cpp =================================================================== --- /issm/trunk-jpl/src/mex/UpdateDynamicConstraints/UpdateDynamicConstraints.cpp (revision 11678) +++ /issm/trunk-jpl/src/mex/UpdateDynamicConstraints/UpdateDynamicConstraints.cpp (revision 11679) @@ -11,5 +11,5 @@ Nodes *nodes = NULL; Parameters *parameters = NULL; - Vec yg = NULL; + Vector* yg = NULL; /*Boot module: */ @@ -32,5 +32,5 @@ /*Free ressources: */ - VecFree(&yg); + delete yg; delete constraints; delete nodes;