Changeset 17653

Timestamp:

04/04/14 16:24:47 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: separating update tau and d

Location:

issm/trunk-jpl/src/c

Files:

• analyses/StressbalanceAnalysis.cpp (modified) (7 diffs)
• analyses/StressbalanceAnalysis.h (modified) (1 diff)
• solutionsequences/solutionsequence_la_theta.cpp (modified) (4 diffs)

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

@@ -4320 +4320 @@
         xDelete<int>(cs_list);
 }/*}}}*/
-void StressbalanceAnalysis::InputUpdateFromSolutionFSXTH(Elements* elements,Parameters* parameters){/*{{{*/
+void StressbalanceAnalysis::InputUpdateFromSolutionFSXTH_d(Elements* elements,Parameters* parameters){/*{{{*/
 
         /*Intermediaries*/
@@ -4330 +4330 @@
         IssmDouble *xyz_list = NULL;
         IssmDouble  Jdet,r;
-        Gauss*      gauss = NULL;
 
         parameters->FindParam(&r,AugmentedLagrangianREnum);
@@ -4391 +4390 @@
                 }
 
-                gauss=element->NewGauss(5);
+                Gauss* gauss=element->NewGauss(5);
                 for(int ig=gauss->begin();ig<gauss->end();ig++){
                         gauss->GaussPoint(ig);
@@ -4522 +4521 @@
                 }
 
+                /*Clean up*/
+                delete gauss;
+                xDelete<IssmDouble>(xyz_list);
+                xDelete<IssmDouble>(tbasis);
+                xDelete<IssmDouble>(Ke);
+                xDelete<IssmDouble>(pe_xx); xDelete<IssmDouble>(d_xx);
+                xDelete<IssmDouble>(pe_yy); xDelete<IssmDouble>(d_yy);
+                xDelete<IssmDouble>(pe_zz); xDelete<IssmDouble>(d_zz);
+                xDelete<IssmDouble>(pe_xy); xDelete<IssmDouble>(d_xy);
+                xDelete<IssmDouble>(pe_xz); xDelete<IssmDouble>(d_xz);
+                xDelete<IssmDouble>(pe_yz); xDelete<IssmDouble>(d_yz);
+        }
+}/*}}}*/
+void StressbalanceAnalysis::InputUpdateFromSolutionFSXTH_tau(Elements* elements,Parameters* parameters){/*{{{*/
+
+        /*Intermediaries*/
+        int         dim,tausize,meshtype;
+        IssmDouble  epsxx,epsyy,epszz,epsxy,epsxz,epsyz,D_scalar;
+        IssmDouble  d_xx,d_yy,d_zz,d_xy,d_xz,d_yz;
+        IssmDouble  sigmapxx,sigmapyy,sigmapzz,sigmapxy,sigmapxz,sigmapyz;
+        IssmDouble  dvx[3],dvy[3],dvz[3];
+        IssmDouble *xyz_list = NULL;
+        IssmDouble  Jdet,r;
+
+        parameters->FindParam(&r,AugmentedLagrangianREnum);
+        parameters->FindParam(&meshtype,MeshTypeEnum);
+        switch(meshtype){
+                case Mesh2DverticalEnum: dim = 2; tausize = 3; break;
+                case Mesh3DEnum:         dim = 3; tausize = 6; break;
+                case Mesh3DtetrasEnum:   dim = 3; tausize = 6; break;
+                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+        }
+
+        for(int i=0;i<elements->Size();i++){
+                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
+
+                /*Get inputs and parameters*/
+                element->GetVerticesCoordinates(&xyz_list);
+                Input* vx_input=element->GetInput(VxEnum);                 _assert_(vx_input);
+                Input* vy_input=element->GetInput(VyEnum);                 _assert_(vy_input);
+                Input* vz_input;
+                if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
+
+                /*Get previous tau*/
+                Input* sigmapxx_input=element->GetInput(DeviatoricStressxxEnum); _assert_(sigmapxx_input);
+                Input* sigmapyy_input=element->GetInput(DeviatoricStressyyEnum); _assert_(sigmapyy_input);
+                Input* sigmapxy_input=element->GetInput(DeviatoricStressxyEnum); _assert_(sigmapxy_input);
+                Input* sigmapzz_input=NULL; Input* sigmapxz_input=NULL; Input* sigmapyz_input=NULL;
+                if(dim==3){
+                        sigmapzz_input=element->GetInput(DeviatoricStresszzEnum); _assert_(sigmapzz_input);
+                        sigmapxz_input=element->GetInput(DeviatoricStressxzEnum); _assert_(sigmapxz_input);
+                        sigmapyz_input=element->GetInput(DeviatoricStressyzEnum); _assert_(sigmapyz_input);
+                }
+
+                /*Get NEW d*/
+                Input* epsxx_input=element->GetInput(StrainRatexxEnum); _assert_(epsxx_input);
+                Input* epsyy_input=element->GetInput(StrainRateyyEnum); _assert_(epsyy_input);
+                Input* epsxy_input=element->GetInput(StrainRatexyEnum); _assert_(epsxy_input);
+                Input* epszz_input=NULL; Input* epsxz_input=NULL; Input* epsyz_input=NULL;
+                if(dim==3){
+                        epszz_input=element->GetInput(StrainRatezzEnum); _assert_(epszz_input);
+                        epsxz_input=element->GetInput(StrainRatexzEnum); _assert_(epsxz_input);
+                        epsyz_input=element->GetInput(StrainRateyzEnum); _assert_(epsyz_input);
+                }
+
+                /*Fetch number of nodes and dof for this finite element*/
+                int tnumnodes = element->GetNumberOfVertices();      //Tensors, P1 DG
+
                 /*Update tau accordingly*/
                 IssmDouble* tau_xx = xNew<IssmDouble>(tnumnodes);
@@ -4529 +4596 @@
                 IssmDouble* tau_xz = NULL;
                 IssmDouble* tau_yz = NULL;
-                delete gauss;
-                gauss = element->NewGauss();
+                Gauss* gauss = element->NewGauss();
                 for(int ig=0;ig<tnumnodes;ig++){
                         gauss->GaussNode(P1DGEnum,ig);
@@ -4559 +4625 @@
                         }
 
+                        /*Get new d*/
+                        epsxx_input->GetInputValue(&d_xx,gauss);
+                        epsyy_input->GetInputValue(&d_yy,gauss);
+                        epsxy_input->GetInputValue(&d_xy,gauss);
+                        if(dim==3){
+                                epszz_input->GetInputValue(&d_zz,gauss);
+                                epsxz_input->GetInputValue(&d_xz,gauss);
+                                epsyz_input->GetInputValue(&d_yz,gauss);
+                        }
+
                         /*Get d and update tau accordingly*/
-                        tau_xx[ig] = sigmapxx + r*(epsxx - d_xx[ig]);
-                        tau_yy[ig] = sigmapyy + r*(epsyy - d_yy[ig]);
-                        tau_xy[ig] = sigmapxy + r*(epsxy - d_xy[ig]);
+                        tau_xx[ig] = sigmapxx + r*(epsxx - d_xx);
+                        tau_yy[ig] = sigmapyy + r*(epsyy - d_yy);
+                        tau_xy[ig] = sigmapxy + r*(epsxy - d_xy);
                         if(dim==3){
-                                tau_zz[ig] = sigmapzz + r*(epszz - d_zz[ig]);
-                                tau_xz[ig] = sigmapxz + r*(epsxz - d_xz[ig]);
-                                tau_yz[ig] = sigmapyz + r*(epsyz - d_yz[ig]);
+                                tau_zz[ig] = sigmapzz + r*(epszz - d_zz);
+                                tau_xz[ig] = sigmapxz + r*(epsxz - d_xz);
+                                tau_yz[ig] = sigmapyz + r*(epsyz - d_yz);
                         }
                 }
@@ -4583 +4659 @@
                 delete gauss;
                 xDelete<IssmDouble>(xyz_list);
-                xDelete<IssmDouble>(tbasis);
-                xDelete<IssmDouble>(Ke);
-                xDelete<IssmDouble>(pe_xx); xDelete<IssmDouble>(d_xx);   xDelete<IssmDouble>(tau_xx);
-                xDelete<IssmDouble>(pe_yy); xDelete<IssmDouble>(d_yy);   xDelete<IssmDouble>(tau_yy);
-                xDelete<IssmDouble>(pe_zz); xDelete<IssmDouble>(d_zz);   xDelete<IssmDouble>(tau_zz);
-                xDelete<IssmDouble>(pe_xy); xDelete<IssmDouble>(d_xy);   xDelete<IssmDouble>(tau_xy);
-                xDelete<IssmDouble>(pe_xz); xDelete<IssmDouble>(d_xz);   xDelete<IssmDouble>(tau_xz);
-                xDelete<IssmDouble>(pe_yz); xDelete<IssmDouble>(d_yz);   xDelete<IssmDouble>(tau_yz);
+                xDelete<IssmDouble>(tau_xx);
+                xDelete<IssmDouble>(tau_yy);
+                xDelete<IssmDouble>(tau_zz);
+                xDelete<IssmDouble>(tau_xy);
+                xDelete<IssmDouble>(tau_xz);
+                xDelete<IssmDouble>(tau_yz);
         }
 }/*}}}*/

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

@@ -81 +81 @@
                 void GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element);
                 void InputUpdateFromSolutionFS(IssmDouble* solution,Element* element);
-                void InputUpdateFromSolutionFSXTH(Elements* elements,Parameters* parameters);
+                void InputUpdateFromSolutionFSXTH_d(Elements* elements,Parameters* parameters);
+                void InputUpdateFromSolutionFSXTH_tau(Elements* elements,Parameters* parameters);
                 void InitializeXTH(Elements* elements,Parameters* parameters);
                 /*Coupling*/

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

@@ -22 +22 @@
         int  configuration_type,max_nonlinear_iterations;
 
+        IssmDouble r     = .6;
+        IssmDouble theta = 0.; // 0<theta<.5   -> .15<theta<.45
+
         /*Create analysis*/
         StressbalanceAnalysis* analysis = new StressbalanceAnalysis();
@@ -30 +33 @@
         femmodel->parameters->FindParam(&max_nonlinear_iterations,StressbalanceMaxiterEnum);
         femmodel->parameters->SetParam(.6,AugmentedLagrangianREnum);
+        femmodel->parameters->SetParam(theta,AugmentedLagrangianThetaEnum);
 
         /*Update constraints and initialize d and tau if necessary*/
@@ -45 +49 @@
                 delete ug_old;ug_old=ug;
 
+                /*Calculate d*/
+                if(theta){
+                        analysis->InputUpdateFromSolutionFSXTH_d(femmodel->elements,femmodel->parameters);
+                }
+
                 /*Solve KU=F*/
                 SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
@@ -57 +66 @@
 
                 /*Update d and tau accordingly*/
-                analysis->InputUpdateFromSolutionFSXTH(femmodel->elements,femmodel->parameters);
+                analysis->InputUpdateFromSolutionFSXTH_d(  femmodel->elements,femmodel->parameters);
+                analysis->InputUpdateFromSolutionFSXTH_tau(femmodel->elements,femmodel->parameters);
 
                 /*Check for convergence*/