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

../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h

@@ -74 +74 @@
                 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);

../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

@@ -944 +944 @@
                 #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 +3290 @@
         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*/