Changeset 7391


Ignore:
Timestamp:
02/08/11 21:10:57 (14 years ago)
Author:
Eric.Larour
Message:

Added support for Stokes solvers

Location:
issm/trunk/src
Files:
2 added
27 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/Makefile.am

    r7312 r7391  
    11391139                                        ./modules/Solverx/Solverx.cpp\
    11401140                                        ./modules/Solverx/Solverx.h\
     1141                                        ./modules/Solverx/DofTypesToIndexSet.cpp\
    11411142                                        ./modules/Scotchx/Scotchx.cpp\
    11421143                                        ./modules/Scotchx/Scotchx.h\

  • issm/trunk/src/c/modules/Solverx/Solverx.cpp

    r7190 r7391  
    1313#endif
    1414
    15 void    Solverx(Vec* puf, Mat Kff, Vec pf, Vec uf0,Parameters* parameters){
     15void    Solverx(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters){
    1616
    1717        /*output: */
    … …  
    2929        int        solver_type;
    3030        bool       fromlocalsize    = true;
     31        /*stokes: */
     32        IS         isv=NULL;
     33        IS         isp=NULL;
    3134
    3235        #if _PETSC_VERSION_ == 3
    … …  
    3538        #endif
    3639
     40
    3741        /*Display message*/
    3842        _printf_(VerboseModule(),"   Solving\n");
    3943        if(VerboseSolver())PetscOptionsPrint(stdout);
    4044
    41         /*First, check that f-set is not NULL, ie model is fully constrained: */
     45        /*First, check that f-set is not NULL, ie model is fully constrained: {{{*/
    4246        _assert_(Kff);
    4347        MatGetSize(Kff,&global_m,&global_n); _assert_(global_m==global_m);
    … …  
    4549                *puf=NULL; return;
    4650        }
    47 
     51        /*}}}*/
     52        /*Initial guess logic here: {{{1*/
    4853        /*Now, check that we are not getting an initial guess to the solver, if we are running a direct solver: */
    4954        #if _PETSC_VERSION_ == 3
    … …  
    6065                MatGetLocalSize(Kff,&local_m,&local_n);uf=NewVec(local_n,fromlocalsize);
    6166        }
    62 
    63         /*Process petsc options to see if we are not using special types of external solvers: */
     67        /*}}}*/
     68        /*Process petsc options to see if we are not using special types of external solvers: {{{1*/
    6469        PetscOptionsDetermineSolverType(&solver_type);
    6570
    … …  
    8893                MatConvert(Kff,MATSUPERLU_DIST,MAT_REUSE_MATRIX,&Kff);
    8994        }
     95        if (solver_type==StokesSolverEnum){
     96                _error_("Petsc 2 does not support multi-physics solvers");
     97        }
    9098        #endif
    9199        #endif
    92 
    93         /*Prepare solver*/
     100        /*}}}*/
     101        /*Prepare solver:{{{1*/
    94102        KSPCreate(MPI_COMM_WORLD,&ksp);
    95103        KSPSetOperators(ksp,Kff,Kff,DIFFERENT_NONZERO_PATTERN);
    96104        KSPSetFromOptions(ksp);
    97105
    98         #ifdef _SERIAL_
    99         #if _PETSC_VERSION_ == 3
     106        #if defined(_SERIAL_) && _PETSC_VERSION_==3
    100107        /*specific solver?: */
    101108        KSPGetPC(ksp,&pc);
    … …  
    104111        }
    105112        #endif
     113
     114        #if defined(_PARALLEL_) && _PETSC_VERSION_==3
     115        /*Stokes: */
     116        if (solver_type==StokesSolverEnum){
     117                /*Make indices out of doftypes: */
     118                if(!df)_error_("need doftypes for Stokes solver!\n");
     119                DofTypesToIndexSet(&isv,&isp,df,StokesSolverEnum);
     120
     121                /*Set field splits: */
     122                KSPGetPC(ksp,&pc);
     123                PCFieldSplitSetIS(pc,isv);
     124                PCFieldSplitSetIS(pc,isp);
     125
     126        }
    106127        #endif
    107128
    108         /*If initial guess for solution, use it, except if we are using the MUMPS direct solver, where any initial
    109          * guess will crash Petsc: */
     129        /*}}}*/
     130        /*If initial guess for solution, use it, except if we are using the MUMPS direct solver, where any initial guess will crash Petsc: {{{1*/
    110131        if (uf0){
    111132                if( (solver_type!=MUMPSPACKAGE_LU) && (solver_type!=MUMPSPACKAGE_CHOL) && (solver_type!=SPOOLESPACKAGE_LU)&& (solver_type!=SPOOLESPACKAGE_CHOL)&& (solver_type!=SUPERLUDISTPACKAGE)){
    … …  
    113134                }
    114135        }
     136        /*}}}*/
    115137       
    116138        if(VerboseSolver())KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);
    … …  
    125147        /*Free ressources:*/
    126148        KSPFree(&ksp);
    127        
     149               
    128150        /*Assign output pointers:*/
    129151        *puf=uf;

  • issm/trunk/src/c/modules/Solverx/Solverx.h

    r5686 r7391  
    99
    1010/* local prototypes: */
    11 void    Solverx( Vec* puf, Mat Kff, Vec pf, Vec uf0,Parameters* parameters);
    12 
     11void    Solverx( Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df,Parameters* parameters);
     12void    DofTypesToIndexSet(IS* pisv, IS* pisp, Vec df,int typeenum);
    1313#endif  /* _SOLVERX_H */
    1414

  • issm/trunk/src/c/modules/SystemMatricesx/SystemMatricesx.cpp

    r6862 r7391  
    99#include "../../EnumDefinitions/EnumDefinitions.h"
    1010
    11 void SystemMatricesx(Mat* pKgg, Mat* pKff, Mat* pKfs, Vec* ppg, Vec* ppf, double* pkmax,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,bool kflag,bool pflag,bool penalty_kflag,bool penalty_pflag){
     11void SystemMatricesx(Mat* pKgg, Mat* pKff, Mat* pKfs, Vec* ppg, Vec* ppf, Vec* pdg, 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){
    1212       
    1313        /*intermediary: */
    … …  
    2626        Vec    pg   = NULL;
    2727        Vec    pf   = NULL;
     28        Vec    dg=NULL;
     29        Vec    df=NULL;
    2830        double kmax = 0;
    2931
    … …  
    5254        if(kflag){
    5355
    54                 if(!buildkff)Kgg=NewMat(gsize,gsize,connectivity,numberofdofspernode);
     56                if(!buildkff){
     57                        Kgg=NewMat(gsize,gsize,connectivity,numberofdofspernode);
     58                        dg=NewVec(gsize);
     59                }
    5560                else{
    5661                        Kff=NewMat(fsize,fsize,connectivity,numberofdofspernode);
    5762                        Kfs=NewMat(fsize,ssize,connectivity,numberofdofspernode);
     63                        df=NewVec(fsize);
    5864                }
    5965
    … …  
    6167                for (i=0;i<elements->Size();i++){
    6268                        element=(Element*)elements->GetObjectByOffset(i);
    63                         element->CreateKMatrix(Kgg,Kff,Kfs);
     69                        element->CreateKMatrix(Kgg,Kff,Kfs,dg,df);
    6470                }
    6571
    … …  
    7076                }
    7177
    72                 /*Assemble matrix and compress matrix to save memory: */
     78                /*Assemble matrix and doftypes and compress matrix to save memory: */
    7379                if(!buildkff){
    7480                        MatAssemblyBegin(Kgg,MAT_FINAL_ASSEMBLY);
    … …  
    7783                        MatCompress(Kgg);
    7884                        #endif
     85                        VecAssemblyBegin(dg);
     86                        VecAssemblyEnd(dg);
    7987                }
    8088                else{
    … …  
    9098                        MatCompress(Kfs);
    9199                        #endif
     100                        VecAssemblyBegin(df);
     101                        VecAssemblyEnd(df);
    92102                }
     103               
    93104        }
    94 if(pflag){
     105       
     106        if(pflag){
    95107
    96108                if(!buildkff)pg=NewVec(gsize);
    … …  
    154166                }
    155167        }
    156 if(penalty_pflag){
     168
     169       
     170        if(penalty_pflag){
    157171
    158172                /*Fill right hand side vector, from loads: */
    … …  
    178192        if(ppg)  *ppg=pg;
    179193        if(ppf)  *ppf=pf;
     194        if(pdg)  *pdg=dg;
     195        if(pdf)  *pdf=df;
    180196        if(pkmax) *pkmax=kmax;
    181197}

  • issm/trunk/src/c/modules/SystemMatricesx/SystemMatricesx.h

    r6010 r7391  
    1010
    1111/* local prototypes: */
    12 void SystemMatricesx(Mat* pKgg, Mat* pKff, Mat* pKfs, Vec* ppg, Vec* ppf, double* pkmax,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,
     12void SystemMatricesx(Mat* pKgg, Mat* pKff, Mat* pKfs, Vec* ppg, Vec* ppf, Vec* pdg, Vec* pdf, double* pkmax,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,
    1313                        bool kflag=true,bool pflag=true,bool penalty_kflag=true,bool penalty_pflag=true);
    1414

  • issm/trunk/src/c/objects/Elements/Element.h

    r7323 r7391  
    2828                virtual void   Configure(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters)=0;
    2929                virtual void   SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters)=0;
    30                 virtual void   CreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs)=0;
     30                virtual void   CreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs,Vec dg, Vec df)=0;
    3131                virtual void   CreatePVector(Vec pg, Vec pf)=0;
    3232                virtual void   GetSolutionFromInputs(Vec solution)=0;

  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r7374 r7391  
    529529}/*}}}*/
    530530/*FUNCTION Penta::CreateKMatrix {{{1*/
    531 void  Penta::CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs){
     531void  Penta::CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs,Vec dg, Vec df){
    532532
    533533        /*retrieve parameters: */
    534534        ElementMatrix* Ke=NULL;
     535        ElementVector* De=NULL;
    535536        int analysis_type;
    536537        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    … …  
    546547        switch(analysis_type){
    547548                case DiagnosticHorizAnalysisEnum: case AdjointHorizAnalysisEnum:
    548                         Ke=CreateKMatrixDiagnosticHoriz();
     549                        Ke=CreateKMatrixDiagnosticHoriz(); De=CreateDVectorDiagnosticHoriz();
    549550                        break;
    550551                case DiagnosticHutterAnalysisEnum:
    … …  
    581582                delete Ke;
    582583        }
     584        /*Add to global Vector*/
     585        if(De){
     586                De->InsertIntoGlobal(dg,df);
     587                delete De;
     588        }
    583589}
    584590/*}}}*/
    … …  
    10061012}
    10071013/*}}}*/
     1014/*FUNCTION Penta::CreateDVectorDiagnosticHoriz {{{1*/
     1015ElementVector* Penta::CreateDVectorDiagnosticHoriz(void){
     1016
     1017        int approximation;
     1018        inputs->GetParameterValue(&approximation,ApproximationEnum);
     1019
     1020        switch(approximation){
     1021                case StokesApproximationEnum:
     1022                        return CreateDVectorDiagnosticStokes();
     1023                default:
     1024                        return NULL; //no need for doftypes outside of stokes approximation
     1025        }
     1026}
     1027/*}}}*/
    10081028/*FUNCTION Penta::CreateKMatrixDiagnosticHutter{{{1*/
    10091029ElementMatrix* Penta::CreateKMatrixDiagnosticHutter(void){
    … …  
    13651385        delete Ke2;
    13661386        return Ke;
     1387}
     1388/*}}}*/
     1389/*FUNCTION Penta::CreateDVectorDiagnosticStokes{{{1*/
     1390ElementVector* Penta::CreateDVectorDiagnosticStokes(void){
     1391
     1392        /*output: */
     1393        ElementVector* De=NULL;
     1394        /*intermediary: */
     1395        int approximation;
     1396        int i;
     1397
     1398        /*Initialize Element vector and return if necessary*/
     1399        inputs->GetParameterValue(&approximation,ApproximationEnum);
     1400        if(approximation!=StokesApproximationEnum) return NULL;
     1401
     1402        De=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     1403
     1404        for (i=0;i<NUMVERTICES;i++){
     1405                De->values[i*4+0]=VelocityEnum;
     1406                De->values[i*4+1]=VelocityEnum;
     1407                De->values[i*4+2]=VelocityEnum;
     1408                De->values[i*4+3]=PressureEnum;
     1409        }
     1410
     1411        return De;
    13671412}
    13681413/*}}}*/

  • issm/trunk/src/c/objects/Elements/Penta.h

    r7323 r7391  
    8080                void   SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters);
    8181                double RegularizeInversion(void);
    82                 void   CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs);
     82                void   CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs,Vec dg,Vec df);
    8383                void   CreatePVector(Vec pg, Vec pf);
    8484                void   DeleteResults(void);
    … …  
    147147                ElementMatrix* CreateKMatrixCouplingPattynStokes(void);
    148148                ElementMatrix* CreateKMatrixDiagnosticHoriz(void);
     149                ElementVector* CreateDVectorDiagnosticHoriz(void);
     150                ElementVector* CreateDVectorDiagnosticStokes(void);
    149151                ElementMatrix* CreateKMatrixDiagnosticHutter(void);
    150152                ElementMatrix* CreateKMatrixDiagnosticMacAyeal2d(void);

  • issm/trunk/src/c/objects/Elements/Tria.cpp

    r7374 r7391  
    377377/*}}}*/
    378378/*FUNCTION Tria::CreateKMatrix {{{1*/
    379 void  Tria::CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs){
     379void  Tria::CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs,Vec dg,Vec df){
    380380
    381381        /*retreive parameters: */

  • issm/trunk/src/c/objects/Elements/Tria.h

    r7323 r7391  
    7777                void   SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters);
    7878                double RegularizeInversion(void);
    79                 void   CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs);
     79                void   CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs,Vec dg,Vec df);
    8080                void   CreatePVector(Vec pg, Vec pf);
    8181                int    GetNodeIndex(Node* node);

  • issm/trunk/src/c/objects/Numerics/ElementVector.cpp

    r6412 r7391  
    197197}
    198198/*}}}*/
     199/*FUNCTION ElementVector::InsertIntoGlobal(Vec pg, Vec pf){{{1*/
     200void ElementVector::InsertIntoGlobal(Vec pg, Vec pf){
     201
     202        int i;
     203        double* localvalues=NULL;
     204
     205        if(!pf){
     206                /*add local values into global  vector, using the fglobaldoflist: */
     207                VecSetValues(pg,this->nrows,this->gglobaldoflist,(const double*)values,INSERT_VALUES);
     208        }
     209        else{
     210                if(this->fsize){
     211                        /*first, retrieve values that are in the f-set from the g-set values vector: */
     212                        localvalues=(double*)xmalloc(this->fsize*sizeof(double));
     213                        for(i=0;i<this->fsize;i++){
     214                                localvalues[i]=this->values[this->flocaldoflist[i]];
     215                        }
     216                        /*add local values into global  vector, using the fglobaldoflist: */
     217                        VecSetValues(pf,this->fsize,this->fglobaldoflist,(const double*)localvalues,INSERT_VALUES);
     218
     219                        /*Free ressources:*/
     220                        xfree((void**)&localvalues);
     221                }
     222        }
     223}
     224/*}}}*/
    199225/*FUNCTION ElementVector::Echo{{{1*/
    200226void ElementVector::Echo(void){

  • issm/trunk/src/c/objects/Numerics/ElementVector.h

    r6027 r7391  
    4242                /*ElementVector specific routines {{{1*/
    4343                void AddToGlobal(Vec pg, Vec pf);
     44                void InsertIntoGlobal(Vec pg, Vec pf);
    4445                void Echo(void);
    4546                void Init(ElementVector* pe);

  • issm/trunk/src/c/solvers/solver_adjoint_linear.cpp

    r6580 r7391  
    1616        Vec  ug  = NULL, uf  = NULL;
    1717        Vec  pg  = NULL, pf  = NULL;
     18        Vec  dg  = NULL, df  = NULL;
    1819        bool kffpartitioning = false;
    1920
    … …  
    2223
    2324        if(kffpartitioning){
    24                 SystemMatricesx(NULL,&Kff, &Kfs, NULL,&pf, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     25                SystemMatricesx(NULL,&Kff, &Kfs, NULL,&pf, NULL, &df, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    2526                Reduceloadx(pf, Kfs, femmodel->ys,true); MatFree(&Kfs); //true means spc = 0
    2627        }
    2728        else{
    28                 SystemMatricesx(&Kgg, NULL, NULL, &pg,NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     29                SystemMatricesx(&Kgg, NULL, NULL, &pg,NULL, &dg, NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    2930                Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets); MatFree(&Kgg);
    3031                Reduceloadfromgtofx(&pf, pg, Kfs, femmodel->ys, femmodel->nodesets,true);VecFree(&pg); MatFree(&Kfs);//true means spc = 0
     32                Reducevectorgtofx(&df, dg, femmodel->nodesets,femmodel->parameters); VecFree(&dg);
    3133        }
    3234
    33         Solverx(&uf, Kff, pf, NULL, femmodel->parameters); MatFree(&Kff); VecFree(&pf);
     35        Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); MatFree(&Kff); VecFree(&pf); VecFree(&df);
    3436        Mergesolutionfromftogx(&ug, uf,femmodel->ys,femmodel->nodesets,true); VecFree(&uf);//true means spc0
    3537        InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); VecFree(&ug); VecFree(&uf);

  • issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp

    r6580 r7391  
    1616        Vec ug  = NULL, uf  = NULL, old_ug= NULL, old_uf = NULL;
    1717        Vec pg  = NULL, pf  = NULL;
     18        Vec dg  = NULL, df  = NULL;
     19       
    1820        Loads* loads=NULL;
    1921        int converged;
    … …  
    5355
    5456                if(kffpartitioning){
    55                         SystemMatricesx(NULL,&Kff, &Kfs, NULL,&pf, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters);
     57                        SystemMatricesx(NULL,&Kff, &Kfs, NULL,&pf, NULL, &df, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters);
    5658                        Reduceloadx(pf, Kfs, femmodel->ys); MatFree(&Kfs);
    5759                }
    5860                else{
    59                         SystemMatricesx(&Kgg, NULL, NULL, &pg,NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters);
     61                        SystemMatricesx(&Kgg, NULL, NULL, &pg,NULL, &dg, NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters);
    6062                        Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets); MatFree(&Kgg);
    6163                        Reduceloadfromgtofx(&pf, pg, Kfs, femmodel->ys, femmodel->nodesets); VecFree(&pg); MatFree(&Kfs);
     64                        Reducevectorgtofx(&df, dg, femmodel->nodesets,femmodel->parameters); VecFree(&dg);
    6265                }
    6366               
    64                 Solverx(&uf, Kff, pf, old_uf, femmodel->parameters);
     67                Solverx(&uf, Kff, pf, old_uf, df, femmodel->parameters);
    6568                Mergesolutionfromftogx(&ug, uf,femmodel->ys,femmodel->nodesets);
    6669                InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug);
    … …  
    6972                _printf_(VerboseConvergence(),"   number of unstable constraints: %i\n",num_unstable_constraints);
    7073
    71                 convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters); MatFree(&Kff);VecFree(&pf);
     74                convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters); MatFree(&Kff);VecFree(&pf); VecFree(&df);
    7275                InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum);
    7376

  • issm/trunk/src/c/solvers/solver_linear.cpp

    r6580 r7391  
    1414        Vec  ug  = NULL, uf  = NULL;
    1515        Vec  pg  = NULL, pf  = NULL;
     16        Vec dg  = NULL, df  = NULL;
    1617        bool kffpartitioning;
    1718
    … …  
    2021
    2122        if(kffpartitioning){
    22                 SystemMatricesx(NULL,&Kff, &Kfs, NULL,&pf, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     23                SystemMatricesx(NULL,&Kff, &Kfs, NULL,&pf, NULL, &df, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    2324                Reduceloadx(pf, Kfs, femmodel->ys); MatFree(&Kfs);
    2425        }
    2526        else{
    26                 SystemMatricesx(&Kgg, NULL, NULL, &pg,NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     27                SystemMatricesx(&Kgg, NULL, NULL, &pg,NULL, &dg, NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    2728                Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets); MatFree(&Kgg);
    2829                Reduceloadfromgtofx(&pf, pg, Kfs, femmodel->ys, femmodel->nodesets);VecFree(&pg); MatFree(&Kfs);
     30                Reducevectorgtofx(&df, dg, femmodel->nodesets,femmodel->parameters); VecFree(&dg);
    2931        }
    3032
    31         Solverx(&uf, Kff, pf, NULL, femmodel->parameters); MatFree(&Kff); VecFree(&pf);
     33        Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); MatFree(&Kff); VecFree(&pf); VecFree(&df);
    3234        Mergesolutionfromftogx(&ug, uf,femmodel->ys,femmodel->nodesets);VecFree(&uf);
    3335        InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); VecFree(&ug); VecFree(&uf);

  • issm/trunk/src/c/solvers/solver_stokescoupling_nonlinear.cpp

    r6580 r7391  
    1616        Vec ug_horiz  = NULL, uf_horiz  = NULL, old_ug_horiz= NULL, old_uf_horiz = NULL;
    1717        Vec pg_horiz  = NULL, pf_horiz  = NULL;
     18        Vec dg_horiz  = NULL, df_horiz  = NULL;
    1819        Mat Kgg_vert  = NULL, Kff_vert  = NULL, Kfs_vert    = NULL;
    1920        Vec ug_vert   = NULL, uf_vert   = NULL;
    2021        Vec pg_vert   = NULL, pf_vert   = NULL;
     22        Vec dg_vert   = NULL, df_vert   = NULL;
    2123        int converged;
    2224        int constraints_converged;
    … …  
    5658
    5759                if(kffpartitioning){
    58                         SystemMatricesx(NULL,&Kff_horiz, &Kfs_horiz, NULL,&pf_horiz, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     60                        SystemMatricesx(NULL,&Kff_horiz, &Kfs_horiz, NULL,&pf_horiz, NULL, &df_horiz, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    5961                        Reduceloadx(pf_horiz, Kfs_horiz, femmodel->ys); MatFree(&Kfs_horiz);
    6062                }
    6163                else{
    62                         SystemMatricesx(&Kgg_horiz, NULL, NULL, &pg_horiz,NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     64                        SystemMatricesx(&Kgg_horiz, NULL, NULL, &pg_horiz,NULL, &dg_horiz,NULL,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    6365                        Reducematrixfromgtofx(&Kff_horiz,&Kfs_horiz,Kgg_horiz,femmodel->nodesets); MatFree(&Kgg_horiz);
    6466                        Reduceloadfromgtofx(&pf_horiz, pg_horiz, Kfs_horiz, femmodel->ys, femmodel->nodesets); VecFree(&pg_horiz); MatFree(&Kfs_horiz);
     67                        Reducevectorgtofx(&df_horiz, dg_horiz, femmodel->nodesets,femmodel->parameters); VecFree(&dg_horiz);
    6568                }
    6669               
    67                 Solverx(&uf_horiz, Kff_horiz, pf_horiz, old_uf_horiz, femmodel->parameters);
     70                Solverx(&uf_horiz, Kff_horiz, pf_horiz, old_uf_horiz, df_horiz,femmodel->parameters);
    6871                Mergesolutionfromftogx(&ug_horiz, uf_horiz,femmodel->ys,femmodel->nodesets);
    6972                InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug_horiz);
    7073
    71                 convergence(&converged,Kff_horiz,pf_horiz,uf_horiz,old_uf_horiz,femmodel->parameters); MatFree(&Kff_horiz);VecFree(&pf_horiz);
     74                convergence(&converged,Kff_horiz,pf_horiz,uf_horiz,old_uf_horiz,femmodel->parameters); MatFree(&Kff_horiz);VecFree(&pf_horiz); VecFree(&df_horiz);
    7275
    7376                /*Second compute vertical velocity: */
    7477                femmodel->SetCurrentConfiguration(DiagnosticVertAnalysisEnum);
    7578                if(kffpartitioning){
    76                         SystemMatricesx(NULL,&Kff_vert, &Kfs_vert, NULL,&pf_vert, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     79                        SystemMatricesx(NULL,&Kff_vert, &Kfs_vert, NULL,&pf_vert, NULL, &df_vert,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    7780                        Reduceloadx(pf_vert, Kfs_vert, femmodel->ys); MatFree(&Kfs_vert);
    7881                }
    7982                else{
    80                         SystemMatricesx(&Kgg_vert, NULL, NULL, &pg_vert,NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     83                        SystemMatricesx(&Kgg_vert, NULL, NULL, &pg_vert,NULL, &dg_horiz, NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    8184                        Reducematrixfromgtofx(&Kff_vert,&Kfs_vert,Kgg_vert,femmodel->nodesets); MatFree(&Kgg_vert);
    8285                        Reduceloadfromgtofx(&pf_vert, pg_vert, Kfs_vert, femmodel->ys, femmodel->nodesets);VecFree(&pg_vert); MatFree(&Kfs_vert);
     86                        Reducevectorgtofx(&df_vert, dg_vert, femmodel->nodesets,femmodel->parameters); VecFree(&dg_vert);
    8387                }
    8488
    85                 Solverx(&uf_vert, Kff_vert, pf_vert, NULL, femmodel->parameters); MatFree(&Kff_vert); VecFree(&pf_vert);
     89                Solverx(&uf_vert, Kff_vert, pf_vert, NULL, df_vert,femmodel->parameters); MatFree(&Kff_vert); VecFree(&pf_vert); VecFree(&df_vert);
    8690                Mergesolutionfromftogx(&ug_vert, uf_vert,femmodel->ys,femmodel->nodesets);VecFree(&uf_vert);
    8791                InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug_vert); VecFree(&ug_vert); VecFree(&uf_vert);

  • issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp

    r6580 r7391  
    2222        Vec pg=NULL;
    2323        Vec pf=NULL;
     24        Vec dg=NULL;
     25        Vec df=NULL;
    2426
    2527        bool converged;
    … …  
    5052
    5153                if(kffpartitioning){
    52                         SystemMatricesx(NULL,&Kff, &Kfs, NULL,&pf,&melting_offset,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     54                        SystemMatricesx(NULL,&Kff, &Kfs, NULL,&pf,NULL, &df, &melting_offset,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    5355                        Reduceloadx(pf, Kfs, femmodel->ys); MatFree(&Kfs);
    5456                }
    5557                else{
    56                         SystemMatricesx(&Kgg,NULL, NULL, &pg,NULL, &melting_offset,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     58                        SystemMatricesx(&Kgg,NULL, NULL, &pg,NULL, &dg, NULL, &melting_offset,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    5759                        Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets); MatFree(&Kgg);
    5860                        Reduceloadfromgtofx(&pf, pg, Kfs, femmodel->ys, femmodel->nodesets);VecFree(&pg); MatFree(&Kfs);
     61                        Reducevectorgtofx(&df, dg, femmodel->nodesets,femmodel->parameters); VecFree(&dg);
    5962                }
    6063
    6164                VecFree(&tf);
    62                 Solverx(&tf, Kff, pf,tf_old, femmodel->parameters);
     65                Solverx(&tf, Kff, pf,tf_old, df, femmodel->parameters);
    6366                VecFree(&tf_old); VecDuplicatePatch(&tf_old,tf);
    64                 MatFree(&Kff);VecFree(&pf);VecFree(&tg);
     67                MatFree(&Kff);VecFree(&pf);VecFree(&tg); VecFree(&df);
    6568
    6669                Mergesolutionfromftogx(&tg, tf,femmodel->ys,femmodel->nodesets);

  • issm/trunk/src/c/toolkits/petsc/patches/PetscOptionsDetermineSolverType.cpp

    r6014 r7391  
    4545                solver_type=SUPERLUDISTPACKAGE;
    4646        }
     47       
     48        PetscOptionsGetString(PETSC_NULL,"-issm_option_solver",&option[0],100,&flag);
     49        if (strcmp(option,"stokes")==0){
     50                solver_type=StokesSolverEnum;
     51        }
    4752
    4853        *psolver_type=solver_type;

  • issm/trunk/src/m/solvers/solver_adjoint_linear.m

    r6581 r7391  
    99
    1010        if kffpartitioning,
    11                 [K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
     11                [K_gg,K_ff,K_fs,p_g,p_f,dg,df,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
    1212                p_f = Reduceload( p_f, K_fs, femmodel.ys,true);
     13
    1314        else
    14                 [K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
     15                [K_gg,K_ff,K_fs,p_g,p_f,dg,df,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
    1516                [K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.nodesets);
    1617                p_f = Reduceloadfromgtof( p_g, K_fs, femmodel.ys, femmodel.nodesets,true);
     18                df=Reducevectorgtof( dg, femmodel.nodesets,femmodel.parameters);
    1719        end
    1820
    1921        issmprintf(VerboseSolver(),'%s%g','      condition number of stiffness matrix: ',condest(K_ff));
    20         u_f=Solver(K_ff,p_f,[],femmodel.parameters);
     22        u_f=Solver(K_ff,p_f,[],df,femmodel.parameters);
    2123        u_g= Mergesolutionfromftog( u_f, femmodel.ys, femmodel.nodesets,true);
    2224

  • issm/trunk/src/m/solvers/solver_diagnostic_nonlinear.m

    r6581 r7391  
    2929
    3030                if kffpartitioning,
    31                         [K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters);
     31                        [K_gg,K_ff,K_fs,p_g,p_f,dg,df,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters);
    3232                        p_f = Reduceload( p_f, K_fs, femmodel.ys);
    3333                else
    34                         [K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters);
     34                        [K_gg,K_ff,K_fs,p_g,p_f,dg,df,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters);
    3535                        [K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.nodesets);
    3636                        p_f = Reduceloadfromgtof( p_g, K_fs, femmodel.ys, femmodel.nodesets);
     37                        df=Reducevectorgtof( dg, femmodel.nodesets,femmodel.parameters);
    3738                end
    3839
    39                 uf=Solver(K_ff,p_f,old_uf,femmodel.parameters);
     40                uf=Solver(K_ff,p_f,old_uf,df,femmodel.parameters);
    4041                ug= Mergesolutionfromftog( uf, femmodel.ys, femmodel.nodesets);
    4142

  • issm/trunk/src/m/solvers/solver_linear.m

    r6581 r7391  
    99
    1010        if kffpartitioning,
    11                 [K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
     11                [K_gg,K_ff,K_fs,p_g,p_f,dg,df,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
    1212                p_f = Reduceload( p_f, K_fs, femmodel.ys);
    1313        else
    14                 [K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
     14                [K_gg,K_ff,K_fs,p_g,p_f,dg,df,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
    1515                [K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.nodesets);
    1616                p_f = Reduceloadfromgtof( p_g, K_fs, femmodel.ys, femmodel.nodesets);
     17                df=Reducevectorgtof( dg, femmodel.nodesets,femmodel.parameters);
    1718        end
    1819
    1920        issmprintf(VerboseSolver(),'%s%g','      condition number of stiffness matrix: ',condest(K_ff));
    20         u_f=Solver(K_ff,p_f,[],femmodel.parameters);
     21        u_f=Solver(K_ff,p_f,[],df,femmodel.parameters);
    2122        u_g= Mergesolutionfromftog( u_f, femmodel.ys, femmodel.nodesets);
    2223

  • issm/trunk/src/m/solvers/solver_stokescoupling_nonlinear.m

    r6589 r7391  
    1 function femmodel=solver_couplingstokes_nonlinear(femmodel,conserve_loads)
     1fpunction femmodel=solver_couplingstokes_nonlinear(femmodel,conserve_loads)
    22%SOLVER_COUPLINGSTOKES_NONLINEAR - core solver of coupling run
    33%
    … …  
    2929
    3030                if kffpartitioning,
    31                         [K_gg_horiz,K_ff_horiz,K_fs_horiz,p_g_horiz,p_f_horiz,kmax_horiz]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
     31                        [K_gg_horiz,K_ff_horiz,K_fs_horiz,p_g_horiz,p_f_horiz,d_g_horiz,d_f_horiz,kmax_horiz]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
    3232                        p_f_horiz = Reduceload( p_f_horiz, K_fs_horiz, femmodel.ys);
    3333                else
    34                         [K_gg_horiz,K_ff_horiz,K_fs_horiz,p_g_horiz,p_f_horiz,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
     34                        [K_gg_horiz,K_ff_horiz,K_fs_horiz,p_g_horiz,p_f_horiz,d_g_horiz,d_f_horiz,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
    3535                        [K_ff_horiz, K_fs_horiz] = Reducematrixfromgtof( K_gg_horiz, femmodel.nodesets);
    3636                        p_f_horiz = Reduceloadfromgtof( p_g_horiz, K_fs_horiz, femmodel.ys, femmodel.nodesets);
     37                        d_f_horiz=Reducevectorgtof( d_g_horiz, femmodel.nodesets,femmodel.parameters);
     38
    3739                end
    3840
    39                 uf_horiz=Solver(K_ff_horiz,p_f_horiz,old_uf_horiz,femmodel.parameters);
     41                uf_horiz=Solver(K_ff_horiz,p_f_horiz,old_uf_horiz,d_f_horiz,femmodel.parameters);
    4042                ug_horiz= Mergesolutionfromftog( uf_horiz, femmodel.ys, femmodel.nodesets);
    4143
    … …  
    4951
    5052                if kffpartitioning,
    51                         [K_gg_vert,K_ff_vert,K_fs_vert,p_g_vert,p_f_vert,kmax_vert]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
     53                        [K_gg_vert,K_ff_vert,K_fs_vert,p_g_vert,p_f_vert,d_g_vert,d_f_vert,kmax_vert]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
    5254                        p_f_vert = Reduceload( p_f_vert, K_fs_vert, femmodel.ys);
    5355                else
    54                         [K_gg_vert,K_ff_vert,K_fs_vert,p_g_vert,p_f_vert,kmax_vert]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
     56                        [K_gg_vert,K_ff_vert,K_fs_vert,p_g_vert,p_f_vert,d_g_vert,d_f_vert,kmax_vert]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
    5557                        [K_ff_vert, K_fs_vert] = Reducematrixfromgtof( K_gg_vert, femmodel.nodesets);
    5658                        p_f_vert = Reduceloadfromgtof( p_g_vert, K_fs_vert, femmodel.ys, femmodel.nodesets);
     59                        d_f_vert=Reducevectorgtof( d_g_vert, femmodel.nodesets,femmodel.parameters);
    5760                end
    5861
    59                 uf_vert=Solver(K_ff_vert,p_f_vert,[],femmodel.parameters);
     62                uf_vert=Solver(K_ff_vert,p_f_vert,[],d_f_vert,femmodel.parameters);
    6063                ug_vert= Mergesolutionfromftog( uf_vert, femmodel.ys, femmodel.nodesets);
    6164

  • issm/trunk/src/m/solvers/solver_thermal_nonlinear.m

    r6581 r7391  
    2424
    2525                if kffpartitioning,
    26                         [K_gg,K_ff,K_fs,p_g,p_f,melting_offset]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
     26                        [K_gg,K_ff,K_fs,p_g,p_f,d_g,d_f,melting_offset]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
    2727                        p_f = Reduceload( p_f, K_fs, femmodel.ys);
    2828                else
    29                         [K_gg,K_ff,K_fs,p_g,p_f,melting_offset]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
     29                        [K_gg,K_ff,K_fs,p_g,p_f,d_g,d_f,melting_offset]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
    3030                        [K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.nodesets);
    3131                        p_f = Reduceloadfromgtof( p_g, K_fs, femmodel.ys, femmodel.nodesets);
     32                        d_f=Reducevectorgtof( d_g, femmodel.nodesets,femmodel.parameters);
     33
    3234                end
    3335
    3436                issmprintf(VerboseSolver(),'%s%g','   condition number of stiffness matrix: ',condest(K_ff));
    35                 t_f=Solver(K_ff,p_f,[],femmodel.parameters);
     37                t_f=Solver(K_ff,p_f,[],d_f,femmodel.parameters);
    3638                t_g= Mergesolutionfromftog( t_f, femmodel.ys, femmodel.nodesets);
    3739

  • issm/trunk/src/mex/Solver/Solver.cpp

    r6412 r7391  
    1212        Vec         uf0           = NULL;
    1313        Vec         uf            = NULL;
     14        Vec         df            = NULL;
    1415        Parameters *parameters    = NULL;
    1516        int         verbose;
    … …  
    4748                FetchData(&pf,PF);
    4849                FetchData(&uf0,UF0);
     50                FetchData(&df,DF);
    4951
    5052                /*Core module: */
    51                 Solverx(&uf, Kff, pf, uf0, parameters);
     53                Solverx(&uf, Kff, pf, uf0, df,parameters);
    5254
    5355                /*Write output*/
    … …  
    6769        VecFree(&uf0);
    6870        VecFree(&uf);
     71        VecFree(&df);
    6972        delete parameters;
    7073
    … …  
    7679{
    7780        _printf_(true,"\n");
    78         _printf_(true,"   usage: [uf] = %s(Kff,pf,uf0,parameters);\n",__FUNCT__);
     81        _printf_(true,"   usage: [uf] = %s(Kff,pf,uf0,df,parameters);\n",__FUNCT__);
    7982        _printf_(true,"\n");
    8083}

  • issm/trunk/src/mex/Solver/Solver.h

    r4236 r7391  
    2323#define PF (mxArray*)prhs[1]
    2424#define UF0 (mxArray*)prhs[2]
    25 #define PARAMETERS (mxArray*)prhs[3]
     25#define DF (mxArray*)prhs[3]
     26#define PARAMETERS (mxArray*)prhs[4]
    2627
    2728/* serial output macros: */
    … …  
    3233#define NLHS  1
    3334#undef NRHS
    34 #define NRHS  4
     35#define NRHS  5
    3536
    3637

  • issm/trunk/src/mex/SystemMatrices/SystemMatrices.cpp

    r6412 r7391  
    2222        Vec    pg   = NULL;
    2323        Vec    pf   = NULL;
     24        Vec    dg   = NULL;
     25        Vec    df   = NULL;
     26
    2427        double kmax;
    2528
    … …  
    5356                FetchData(&penalty_kflag,PENALTYKFLAG);
    5457                FetchData(&penalty_pflag,PENALTYPFLAG);
    55                 SystemMatricesx(&Kgg,&Kff,&Kfs,&pg,&pf,&kmax,elements,nodes,vertices,loads,materials,parameters,kflag,pflag,penalty_kflag,penalty_pflag);
     58                SystemMatricesx(&Kgg,&Kff,&Kfs,&pg,&pf,&dg,&df,&kmax,elements,nodes,vertices,loads,materials,parameters,kflag,pflag,penalty_kflag,penalty_pflag);
    5659        }
    5760        else
    58          SystemMatricesx(&Kgg,&Kff,&Kfs,&pg,&pf,&kmax,elements,nodes,vertices,loads,materials,parameters);
     61         SystemMatricesx(&Kgg,&Kff,&Kfs,&pg,&pf,&dg,&df,&kmax,elements,nodes,vertices,loads,materials,parameters);
    5962
    6063        /*write output datasets: */
    … …  
    6467        WriteData(PG,pg);
    6568        WriteData(PF,pf);
     69        WriteData(DG,dg);
     70        WriteData(DF,df);
    6671        WriteData(KMAX,kmax);
    6772       
    … …  
    7883        VecFree(&pg);
    7984        VecFree(&pf);
     85        VecFree(&dg);
     86        VecFree(&df);
    8087
    8188        /*end module: */
    … …  
    8693{
    8794        _printf_(true,"\n");
    88         _printf_(true,"   usage: [Kgg,Kff,Kfs,pg,pf,kmax] = %s(elements,nodes,vertices,loads,materials,parameters);\n",__FUNCT__);
    89         _printf_(true,"   usage: [Kgg,Kff,Kfs,pg,pf,kmax] = %s(elements,nodes,vertices,loads,materials,parameters,kflag,pflag,penalty_kflag,penalty_pflag);\n",__FUNCT__);
     95        _printf_(true,"   usage: [Kgg,Kff,Kfs,pg,pf,dg,df,kmax] = %s(elements,nodes,vertices,loads,materials,parameters);\n",__FUNCT__);
     96        _printf_(true,"   usage: [Kgg,Kff,Kfs,pg,pf,dg,df,kmax] = %s(elements,nodes,vertices,loads,materials,parameters,kflag,pflag,penalty_kflag,penalty_pflag);\n",__FUNCT__);
    9097        _printf_(true,"\n");
    9198}

  • issm/trunk/src/mex/SystemMatrices/SystemMatrices.h

    r6014 r7391  
    3535#define PG   (mxArray**)&plhs[3]
    3636#define PF   (mxArray**)&plhs[4]
    37 #define KMAX (mxArray**)&plhs[5]
     37#define DG   (mxArray**)&plhs[5]
     38#define DF   (mxArray**)&plhs[6]
     39#define KMAX (mxArray**)&plhs[7]
    3840
    3941/* serial arg counts: */
    4042#undef NLHS
    41 #define NLHS  6
     43#define NLHS  8
    4244#undef NRHS
    4345#define NRHS  10
Note: See TracChangeset for help on using the changeset viewer.