source: issm/oecreview/Archive/22819-23185/ISSM-23161-23162.diff

Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 23161)
+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 23162)
@@ -73,6 +73,7 @@
                ElementMatrix* CreateKMatrixFSViscous(Element* element);
                ElementMatrix* CreateKMatrixFSViscousLA(Element* element);
                ElementMatrix* CreateKMatrixFSViscousXTH(Element* element);
+               ElementMatrix* CreatePressureMassMatrix(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 23161)
+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp       (revision 23162)
@@ -3242,6 +3242,56 @@
        xDelete<int>(cs_list);
        return Ke;
 }/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreatePressureMassMatrix(Element* element){/*{{{*/
+
+       /*Intermediaries*/
+       int         i,dim;
+       IssmDouble  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);
+
+       /* 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);
+
+               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*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*/
Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp
===================================================================
--- ../trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp   (revision 23161)
+++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp   (revision 23162)
@@ -10,7 +10,7 @@
 #include "../analyses/analyses.h"
 
 #ifdef _HAVE_PETSC_
-void SchurCGSolver(Vector<IssmDouble>** puf,Mat Kff, Vec pf, Vec uf0,Vec df,Parameters* parameters){/*{{{*/
+void SchurCGSolver(Vector<IssmDouble>** puf,Mat Kff,Mat Mff,Vec pf, Vec uf0,Vec df,Parameters* parameters){/*{{{*/
 
        /*Initialize output*/
        Vec uf = NULL;
@@ -34,6 +34,7 @@
        Vector<IssmDouble>* pf0 = NULL;
        Vector<IssmDouble>* df  = NULL;
        Vector<IssmDouble>* ys  = NULL;
+       Matrix<IssmDouble>* Mff = NULL;
 
        /*parameters:*/
        int max_nonlinear_iterations;
@@ -71,17 +72,32 @@
                CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
                Reduceloadx(pf, Kfs, ys); delete Kfs;
 
+               /*Create mass matrix*/
+               int fsize; Kff->GetSize(&fsize,&fsize);
+               Mff=new Matrix<IssmDouble>(fsize,fsize,100,4);
+               StressbalanceAnalysis* analysis = new StressbalanceAnalysis();
+               /*Get complete stiffness matrix without penalties*/
+               for(int i=0;i<femmodel->elements->Size();i++){
+                       Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                       ElementMatrix* Me = analysis->CreatePressureMassMatrix(element);
+                       if(Me) Me->AddToGlobal(Mff,NULL);
+                       delete Me;
+               }
+               Mff->Assemble();
+               delete analysis;
+
                /*Solve*/
                femmodel->profiler->Start(SOLVER);
                _assert_(Kff->type==PetscMatType); 
                SchurCGSolver(&uf,
                                        Kff->pmatrix->matrix,
+                                       Mff->pmatrix->matrix,
                                        pf->pvector->vector,
                                        old_uf->pvector->vector,
                                        df->pvector->vector,
                                        femmodel->parameters);
                femmodel->profiler->Stop(SOLVER);
-               delete pf0;
+               delete pf0; delete Mff;
 
                /*Merge solution from f set to g set*/
                Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete ys;