source: issm/oecreview/Archive/18296-19100/ISSM-18349-18350.diff

Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp
===================================================================
--- ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp       (revision 18349)
+++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp       (revision 18350)
@@ -111,6 +111,53 @@
        *pdmax = dmax;
        *pLHS  = LHS;
 }/*}}}*/
+void CreateRHS(Vec* pRHS,Mat K,Mat D,Vec Ml,Vec u,IssmDouble theta,IssmDouble deltat,IssmDouble dmax,FemModel* femmodel,int configuration_type){/*{{{*/
+       /*Create Left Hand side of Lower order solution
+        *
+        * RHS = [ML + (1 − theta) deltaT L^n] u^n
+        *
+        * where L = K + D
+        *
+        */
+
+       /*Intermediaries*/
+       Vec         Ku  = NULL;
+       Vec         Du  = NULL;
+       Vec         RHS = NULL;
+       int         dof;
+       IssmDouble  d;
+
+       /*Initialize vectors*/
+       VecDuplicate(u,&Ku);
+       VecDuplicate(u,&Du);
+       VecDuplicate(u,&RHS);
+
+       /*Create RHS = M*u + (1-theta)*deltat*K*u + (1-theta)*deltat*D*u*/
+       MatMult(K,u,Ku);
+       MatMult(D,u,Du);
+       VecPointwiseMult(RHS,Ml,u);
+       VecAXPBYPCZ(RHS,(1-theta)*deltat,(1-theta)*deltat,1,Ku,Du);
+       VecFree(&Ku);
+       VecFree(&Du);
+
+       /*Penalize Dirichlet boundary*/
+       for(int i=0;i<femmodel->constraints->Size();i++){
+               Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
+               if(constraint->InAnalysis(configuration_type)){
+                       constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
+                       d = d*dmax;
+                       if(dof!=-1){
+                               VecSetValues(RHS,1,&dof,(const double*)&d,INSERT_VALUES);
+                       }
+               }
+       }
+       VecAssemblyBegin(RHS);
+       VecAssemblyEnd(  RHS);
+
+       /*Assign output pointer*/
+       *pRHS = RHS;
+
+}/*}}}*/
 #endif
 void solutionsequence_fct(FemModel* femmodel){
 
@@ -120,18 +167,15 @@
        Matrix<IssmDouble>*  Mc = NULL;
        Vector<IssmDouble>*  ug = NULL;
        Vector<IssmDouble>*  uf = NULL;
-       Vector<IssmDouble>*  ys = NULL;
 
-       IssmDouble theta,deltat,dmax=0.;
-       int        dof,ncols,ncols2,ncols3,rstart,rend;
-       int        configuration_type,analysis_type;
-       double     d,mi;
+       IssmDouble theta,deltat,dmax;
+       int        dof,ncols,ncols2,rstart,rend;
+       int        configuration_type;
+       double     d;
        int*       cols  = NULL;
        int*       cols2 = NULL;
-       int*       cols3 = NULL;
        double*    vals  = NULL;
        double*    vals2 = NULL;
-       double*    vals3 = NULL;
 
        /*Create analysis*/
        MasstransportAnalysis* analysis = new MasstransportAnalysis();
@@ -139,7 +183,6 @@
        /*Recover parameters: */
        femmodel->parameters->FindParam(&deltat,TimesteppingTimeStepEnum);
        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
-       femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
        femmodel->UpdateConstraintsx();
        theta = 0.5;
 
@@ -151,6 +194,8 @@
        /*Convert matrices to PETSc matrices*/
        Mat D_petsc  = NULL;
        Mat LHS      = NULL;
+       Vec RHS      = NULL;
+       Vec u        = NULL;
        Mat K_petsc  = K->pmatrix->matrix;
        Vec Ml_petsc = Ml->pvector->vector;
        Mat Mc_petsc = Mc->pmatrix->matrix;
@@ -166,35 +211,9 @@
        GetSolutionFromInputsx(&ug,femmodel);
        Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters);
        delete ug;
-       Vec u  = uf->pvector->vector;
-       Vec Ku = NULL;
-       Vec Du = NULL;
-       Vec RHS = NULL;
-       VecDuplicate(u,&Ku);
-       VecDuplicate(u,&Du);
-       VecDuplicate(u,&RHS);
-       MatMult(K_petsc,u,Ku);
-       MatMult(D_petsc,u,Du);
-       VecPointwiseMult(RHS,Ml_petsc,u);
-       VecAXPBYPCZ(RHS,(1-theta)*deltat,(1-theta)*deltat,1,Ku,Du);// RHS = M*u + (1-theta)*deltat*K*u + (1-theta)*deltat*D*u
-       VecFree(&Ku);
-       VecFree(&Du);
+       CreateRHS(&RHS,K_petsc,D_petsc,Ml_petsc,uf->pvector->vector,theta,deltat,dmax,femmodel,configuration_type);
        delete uf;
 
-       /*Penalize Dirichlet boundary*/
-       for(int i=0;i<femmodel->constraints->Size();i++){
-               Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
-               if(constraint->InAnalysis(analysis_type)){
-                       constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
-                       d = d*dmax;
-                       if(dof!=-1){
-                               VecSetValues(RHS,1,&dof,(const double*)&d,INSERT_VALUES);
-                       }
-               }
-       }
-       VecAssemblyBegin(RHS);
-       VecAssemblyEnd(  RHS);
-
        /*Go solve!*/
        SolverxPetsc(&u,LHS,RHS,NULL,NULL, femmodel->parameters); 
        MatFree(&LHS);
@@ -324,9 +343,8 @@
        uf =new Vector<IssmDouble>(u);
        VecFree(&u);
 
-       Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys;
-       InputUpdateFromSolutionx(femmodel,ug); 
-       delete ug;  
+       InputUpdateFromSolutionx(femmodel,uf); 
+       delete uf;
 
        #else
        _error_("PETSc needs to be installed");