source: issm/oecreview/Archive/23185-23389/ISSM-23198-23199.diff@ 23390

Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 23198)
+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 23199)
@@ -74,6 +74,7 @@
                ElementMatrix* CreateKMatrixFSViscousLA(Element* element);
                ElementMatrix* CreateKMatrixFSViscousXTH(Element* element);
                ElementMatrix* CreatePressureMassMatrix(Element* element);
+               ElementMatrix* CreateSchurPrecondMatrix(Element* element);
                ElementVector* CreatePVectorFS(Element* element);
                ElementVector* CreatePVectorFSFriction(Element* element);
                ElementVector* CreatePVectorFSFront(Element* element);
Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp       (revision 23198)
+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp       (revision 23199)
@@ -944,9 +944,11 @@
                #ifdef _HAVE_PETSC_
                int solver_type;
                PetscOptionsDetermineSolverType(&solver_type);
+
                if(solver_type==FSSolverEnum) is_schur_cg_solver = true;
                #endif
 
+               
                if(is_schur_cg_solver)
                 solutionsequence_schurcg(femmodel);
                else if (fe_FS==XTaylorHoodEnum)
@@ -3288,6 +3290,64 @@
        xDelete<IssmDouble>(pbasis);
        return Ke;
 }/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreateSchurPrecondMatrix(Element* element){/*{{{*/
+
+       /*Intermediaries*/
+       int         i,dim;
+       IssmDouble  viscosity,FSreconditioning,Jdet;
+       IssmDouble *xyz_list = NULL;
+
+       /*Get problem dimension*/
+       element->FindParam(&dim,DomainDimensionEnum);
+
+       /*Fetch number of nodes and dof for this finite element*/
+       int vnumnodes = element->NumberofNodesVelocity();
+       int pnumnodes = element->NumberofNodesPressure();
+       int numdof    = vnumnodes*dim + pnumnodes;
+
+       /*Initialize Element matrix and vectors*/
+       ElementMatrix* Ke  = element->NewElementMatrix(FSvelocityEnum);
+       IssmDouble* pbasis = xNew<IssmDouble>(pnumnodes);
+
+       /*Retrieve all inputs and parameters*/
+       element->GetVerticesCoordinates(&xyz_list);
+       //element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
+       Input* vx_input=element->GetInput(VxEnum);     _assert_(vx_input);
+       Input* vy_input=element->GetInput(VyEnum);     _assert_(vy_input);
+       Input* vz_input = NULL;
+       if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
+
+
+       /* Start  looping on the number of gaussian points: */
+       Gauss* gauss = element->NewGauss(5);
+       for(int ig=gauss->begin();ig<gauss->end();ig++){
+               gauss->GaussPoint(ig);
+
+               element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+               element->NodalFunctionsPressure(pbasis,gauss);
+               element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
+
+
+               if(dim==3){
+                       /*Pressure mass matrix*/
+                       for(int k=0;k<pnumnodes;k++){
+                               for(int j=0;j<pnumnodes;j++){
+                                       Ke->values[(3*vnumnodes+k)*numdof+3*vnumnodes+j] += gauss->weight*1./viscosity*Jdet*(pbasis[j]*pbasis[k]);
+                               }
+                       }
+               }
+               else{
+                       _error_("STOP");
+               }
+       }
+
+       /*Clean up and return*/
+       delete gauss;
+       xDelete<IssmDouble>(xyz_list);
+       xDelete<IssmDouble>(pbasis);
+       return Ke;
+}/*}}}*/
+
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSViscousLA(Element* element){/*{{{*/
 
        /*Intermediaries*/