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

../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp

@@ -111 +111 @@
         *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 +167 @@
         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 +183 @@
         /*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 +194 @@
         /*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 +211 @@
         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 +343 @@
         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");