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

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

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

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

@@ -3242 +3242 @@
         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*/

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

@@ -10 +10 @@
 #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 +34 @@
         Vector<IssmDouble>* pf0 = NULL;
         Vector<IssmDouble>* df  = NULL;
         Vector<IssmDouble>* ys  = NULL;
+        Matrix<IssmDouble>* Mff = NULL;
 
         /*parameters:*/
         int max_nonlinear_iterations;
@@ -71 +72 @@
                 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;