Changeset 18654

Timestamp:

10/17/14 11:57:36 (10 years ago)

Author:

seroussi

Message:

BUG: fixed augmentation for non linear case

Location:

issm/trunk-jpl/src/c

Files:

• analyses/UzawaPressureAnalysis.cpp (modified) (3 diffs)
• solutionsequences/solutionsequence_la.cpp (modified) (4 diffs)

issm/trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp

@@ -81 +81 @@
         IssmDouble*    M    = xNew<IssmDouble>(numnodes);
 
+        IssmDouble connectivity;
         /*Retrieve all inputs and parameters*/
         element->GetVerticesCoordinates(&xyz_list);
-
-        Gauss* gauss = element->NewGauss(5);
-        for(int ig=gauss->begin();ig<gauss->end();ig++){
-                gauss->GaussPoint(ig);
-
-                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-                this->GetM(M,element,gauss);
-                D_scalar=gauss->weight*Jdet;
-                TripleMultiply(M,1,numnodes,1,
-                                        &D_scalar,1,1,0,
-                                        M,1,numnodes,0,
-                                        &Ke->values[0],1);
+        int numvertices = element->GetNumberOfVertices();
+
+        //Gauss* gauss = element->NewGauss(5);
+        //for(int ig=gauss->begin();ig<gauss->end();ig++){
+        //      gauss->GaussPoint(ig);
+
+        //      element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+        //      this->GetM(M,element,gauss);
+        //      D_scalar=gauss->weight*Jdet;
+        //      //TripleMultiply(M,1,numnodes,1,
+        //      //                      &D_scalar,1,1,0,
+        //      //                      M,1,numnodes,0,
+        //      //                      &Ke->values[0],1);
+
+        //}
+        for(int iv=0;iv<numvertices;iv++){
+                connectivity=(IssmDouble)element->VertexConnectivity(iv);
+                Ke->values[iv*numvertices+iv]=1./connectivity;
         }
 
         /*Clean up and return*/
-        delete gauss;
+        //delete gauss;
         xDelete<IssmDouble>(xyz_list);
         xDelete<IssmDouble>(M);
@@ -202 +209 @@
         IssmDouble* valueslambda  = xNewZeroInit<IssmDouble>(numnodessigma);
         IssmDouble* pressure      = xNew<IssmDouble>(numnodes);
-        Input* sigmann_input      = element->GetInput(SigmaNNEnum); _assert_(sigmann_input);
         Input* vx_input           = element->GetInput(VxEnum);      _assert_(vx_input);
         Input* vy_input           = element->GetInput(VyEnum);      _assert_(vy_input);
@@ -217 +223 @@
         /*Now compute sigmann if on base*/
         if(element->IsOnBase() && 0){ 
+                Input* sigmann_input      = element->GetInput(SigmaNNEnum); _assert_(sigmann_input);
                 if(dim==3) _error_("not implemented yet");
 

issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la.cpp

@@ -39 +39 @@
 
         /*Convergence criterion*/
-        int  count = 0;
-        Vector<IssmDouble>* vel     = NULL;
-        Vector<IssmDouble>* vel_old = NULL;
+        int  count       = 0;
+        int  count_local = 0;
+        Vector<IssmDouble>* vel           = NULL;
+        Vector<IssmDouble>* vel_old       = NULL;
+        Vector<IssmDouble>* vel_old_local = NULL;
         GetVectorFromInputsx(&vel,femmodel,VxEnum,VertexEnum);
         GetVectorFromInputsx(&pug,femmodel,PressureEnum,VertexEnum);
@@ -49 +51 @@
 
                 /*save pointer to old velocity*/
-                delete vel_old;vel_old=vel;
+                delete vel_old;vel_old=vel->Duplicate(); vel->Copy(vel_old);
                 delete pug_old;pug_old=pug; 
                 pressure_converged=false; vel_converged=false;
 
-                /*Solve KU=F*/
-                femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
-                femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
-                SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
-                CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
-                Reduceloadx(pf, Kfs, ys); delete Kfs;
-                Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); 
-                delete Kff; delete pf; delete df;
-                Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys;
+                while(true){
+                        count_local++;
+                        delete vel_old_local;vel_old_local=vel;
+                        /*Solve KU=F*/
+                        femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
+                        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+                        SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
+                        CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
+                        Reduceloadx(pf, Kfs, ys); delete Kfs;
+                        Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); 
+                        delete Kff; delete pf; delete df;
+                        Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys;
 
-                /*Update solution*/
-                InputUpdateFromSolutionx(femmodel,ug); delete ug;
-                GetVectorFromInputsx(&vel,femmodel,VxEnum,VertexEnum);
+                        /*Update solution*/
+                        InputUpdateFromSolutionx(femmodel,ug); delete ug;
+                        GetVectorFromInputsx(&vel,femmodel,VxEnum,VertexEnum);
+                        /*Check for convergence*/
+                        Vector<IssmDouble>* dvel_local=vel_old_local->Duplicate(); vel_old_local->Copy(dvel_local); dvel_local->AYPX(vel,-1.0);
+                        IssmDouble ndu=dvel_local->Norm(NORM_TWO);   delete dvel_local;
+                        IssmDouble nu =vel_old_local->Norm(NORM_TWO);
+                        if (xIsNan<IssmDouble>(ndu) || xIsNan<IssmDouble>(nu)) _error_("convergence criterion is NaN!");
+                        if((ndu/nu)<eps_rel){
+                                if(VerboseConvergence()) _printf0_(setw(50) << left << "   Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " < " << eps_rel*100 << " %\n");
+                                break;
+                        }
+                        else{ 
+                                if(VerboseConvergence()) _printf0_(setw(50) << left << "   Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " > " << eps_rel*100 << " %\n");
+                        }
+                        if(count_local>=max_nonlinear_iterations){
+                                _printf0_("   maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n"); 
+                                break;
+                        }
+                }
+                count_local=0;
 
                 femmodel->SetCurrentConfiguration(UzawaPressureAnalysisEnum);
@@ -81 +104 @@
 
                 /*Check for convergence*/
-                Vector<IssmDouble>* dvel=vel_old->Duplicate(); vel_old->Copy(dvel); dvel->AYPX(vel,-1.0);
+                Vector<IssmDouble>* dvel=vel_old_local->Duplicate(); vel_old_local->Copy(dvel); dvel->AYPX(vel,-1.0);
                 IssmDouble ndu=dvel->Norm(NORM_TWO);   delete dvel;
-                IssmDouble nu =vel_old->Norm(NORM_TWO);
+                IssmDouble nu =vel_old_local->Norm(NORM_TWO);
                 if (xIsNan<IssmDouble>(ndu) || xIsNan<IssmDouble>(nu)) _error_("convergence criterion is NaN!");
                 if((ndu/nu)<eps_rel){
@@ -119 +142 @@
         delete vel;  
         delete vel_old;  
+        delete vel_old_local;  
         delete stressanalysis;
         delete pressureanalysis;