Changeset 17540

Timestamp:

03/25/14 13:42:51 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: moved initialization to solution sequence

Location:

issm/trunk-jpl/src/c

Files:

• analyses/StressbalanceAnalysis.cpp (modified) (3 diffs)
• analyses/StressbalanceAnalysis.h (modified) (2 diffs)
• classes/Elements/TriaRef.cpp (modified) (3 diffs)
• solutionsequences/solutionsequence_la_theta.cpp (modified) (2 diffs)

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

@@ -234 +234 @@
                 iomodel->FetchDataToInput(elements,PressureEnum,0.);
                 if(dakota_analysis)elements->InputDuplicate(PressureEnum,QmuPressureEnum);
-                if(finiteelement==XTaylorHoodEnum) this->InitializeXTH(elements);
         }
         if(islevelset){
@@ -821 +820 @@
         }
 }/*}}}*/
-void StressbalanceAnalysis::InitializeXTH(Elements* elements){/*{{{*/
-
-        /*Intermediaries*/
-        int        meshtype,dim;
-        IssmDouble dvx[3],dvy[3],dvz[3];
-        IssmDouble viscosity;
-        IssmDouble *xyz_list = NULL;
-
-        for(int i=0;i<elements->Size();i++){
-                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
-
-                /*Get problem dimension*/
-                element->FindParam(&meshtype,MeshTypeEnum);
-                switch(meshtype){
-                        case Mesh2DverticalEnum: dim = 2; break;
-                        case Mesh3DEnum:         dim = 3; break;
-                        case Mesh3DtetrasEnum:   dim = 3; break;
-                        default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
-                }
-
-                int tnumnodes = element->GetNumberOfVertices();      //Tensors, P1 DG
-
-                /*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 d and tau*/
-                Gauss* gauss = element->NewGauss();
-                for(int i=0;i<tnumnodes;i++){
-                        gauss->GaussNode(P1DGEnum,i);
-
-                        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
-                        vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
-                        if(dim==3){
-                                vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
-                        }
-
-                        _error_("STOP");
-                        element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
-                }
-
-                /*Clean up*/
-                xDelete<IssmDouble>(xyz_list);
-        }
-
-}/*}}}*/
 
 /*Finite Element Analysis*/
@@ -4109 +4059 @@
         xDelete<IssmDouble>(vvalues);
 }/*}}}*/
+void StressbalanceAnalysis::InitializeXTH(Elements* elements,Parameters* parameters){/*{{{*/
+
+        /*Intermediaries*/
+        int        meshtype,dim;
+        IssmDouble dvx[3],dvy[3],dvz[3];
+        IssmDouble viscosity;
+        IssmDouble *xyz_list = NULL;
+
+        /*Get problem dimension*/
+        parameters->FindParam(&meshtype,MeshTypeEnum);
+        switch(meshtype){
+                case Mesh2DverticalEnum: dim = 2; break;
+                case Mesh3DEnum:         dim = 3; break;
+                case Mesh3DtetrasEnum:   dim = 3; 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);}
+
+                /*Allocate new inputs*/
+                int tnumnodes = element->GetNumberOfVertices();      //Tensors, P1 DG
+                IssmDouble* epsxx = xNew<IssmDouble>(tnumnodes); IssmDouble* sigmapxx = xNew<IssmDouble>(tnumnodes);
+                IssmDouble* epsyy = xNew<IssmDouble>(tnumnodes); IssmDouble* sigmapyy = xNew<IssmDouble>(tnumnodes); 
+                IssmDouble* epsxy = xNew<IssmDouble>(tnumnodes); IssmDouble* sigmapxy = xNew<IssmDouble>(tnumnodes); 
+                IssmDouble* epszz = NULL;                        IssmDouble* sigmapzz = NULL;  
+                IssmDouble* epsxz = NULL;                        IssmDouble* sigmapxz = NULL;  
+                IssmDouble* epsyz = NULL;                        IssmDouble* sigmapyz = NULL;  
+                if(dim==3){
+                        epszz = xNew<IssmDouble>(tnumnodes); sigmapzz = xNew<IssmDouble>(tnumnodes);
+                        epsxz = xNew<IssmDouble>(tnumnodes); sigmapxz = xNew<IssmDouble>(tnumnodes); 
+                        epsyz = xNew<IssmDouble>(tnumnodes); sigmapyz = xNew<IssmDouble>(tnumnodes); 
+                }
+
+                /*Get d and tau*/
+                Gauss* gauss = element->NewGauss();
+                for(int i=0;i<tnumnodes;i++){
+                        gauss->GaussNode(P1DGEnum,i);
+
+                        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+                        vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+                        if(dim==3){
+                                vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
+                        }
+
+                        element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
+                        epsxx[i] = dvx[0];                sigmapxx[i] = 2.*viscosity*epsxx[i];
+                        epsyy[i] = dvy[1];                sigmapyy[i] = 2.*viscosity*epsyy[i];
+                        epsxy[i] = 0.5*(dvx[1] + dvy[0]); sigmapxy[i] = 2.*viscosity*epsxy[i];
+                        if(dim==3){
+                                epszz[i] = dvz[2];                sigmapzz[i] = 2.*viscosity*epszz[i];
+                                epsxz[i] = 0.5*(dvx[2] + dvz[0]); sigmapxz[i] = 2.*viscosity*epsxz[i];
+                                epsyz[i] = 0.5*(dvy[2] + dvz[1]); sigmapyz[i] = 2.*viscosity*epsyz[i];
+                        }
+                }
+
+                /*Add inputs*/
+                element->AddInput(StrainRatexxEnum,epsxx,P1DGEnum); element->AddInput(DeviatoricStressxxEnum,sigmapxx,P1DGEnum);
+                element->AddInput(StrainRateyyEnum,epsyy,P1DGEnum); element->AddInput(DeviatoricStressyyEnum,sigmapyy,P1DGEnum);
+                element->AddInput(StrainRatexyEnum,epsxy,P1DGEnum); element->AddInput(DeviatoricStressxyEnum,sigmapxy,P1DGEnum);
+                if(dim==3){
+                        element->AddInput(StrainRatezzEnum,epszz,P1DGEnum); element->AddInput(DeviatoricStresszzEnum,sigmapzz,P1DGEnum);
+                        element->AddInput(StrainRatexzEnum,epsxz,P1DGEnum); element->AddInput(DeviatoricStressxzEnum,sigmapxz,P1DGEnum);
+                        element->AddInput(StrainRateyzEnum,epsyz,P1DGEnum); element->AddInput(DeviatoricStressyzEnum,sigmapyz,P1DGEnum);
+                }
+
+                /*Clean up*/
+                xDelete<IssmDouble>(xyz_list);
+                xDelete<IssmDouble>(epsxx);
+                xDelete<IssmDouble>(epsyy);
+                xDelete<IssmDouble>(epszz);
+                xDelete<IssmDouble>(epsxy);
+                xDelete<IssmDouble>(epsxz);
+                xDelete<IssmDouble>(epsyz);
+        }
+
+}/*}}}*/
 void StressbalanceAnalysis::InputUpdateFromSolutionFS(IssmDouble* solution,Element* element){/*{{{*/
 

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

@@ -19 +19 @@
                 void CreateConstraints(Constraints* constraints,IoModel* iomodel);
                 void CreateLoads(Loads* loads, IoModel* iomodel);
-                void InitializeXTH(Elements* elements);
 
                 /*Finite element Analysis*/
@@ -82 +81 @@
                 void GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element);
                 void InputUpdateFromSolutionFS(IssmDouble* solution,Element* element);
+                void InitializeXTH(Elements* elements,Parameters* parameters);
                 /*Coupling*/
                 ElementMatrix* CreateKMatrixSSA3d(Element* element);

issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp

@@ -246 +246 @@
                         this->element_type = TaylorHoodEnum;
                         return;
+                case XTaylorHoodEnum:
+                        this->element_type = P2Enum;
+                        this->GetNodalFunctions(basis,gauss);
+                        this->element_type = XTaylorHoodEnum;
+                        return;
                 default:
                         _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
@@ -280 +285 @@
                         this->GetNodalFunctions(basis,gauss);
                         this->element_type = TaylorHoodEnum;
+                        return;
+                case XTaylorHoodEnum:
+                        this->element_type = P1Enum;
+                        this->GetNodalFunctions(basis,gauss);
+                        this->element_type = XTaylorHoodEnum;
                         return;
                 default:
@@ -439 +449 @@
                         this->element_type = P2Enum;
                         this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-                        this->element_type = TaylorHoodEnum;
+                        this->element_type = XTaylorHoodEnum;
                         return;
                 default:

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

@@ -20 +20 @@
         int  configuration_type;
 
+        /*Create analysis*/
+        StressbalanceAnalysis* analysis = new StressbalanceAnalysis();
+
         /*Recover parameters: */
         femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+
+        /*Update constraints and initialize d and tau if necessary*/
         femmodel->UpdateConstraintsx();
+        analysis->InitializeXTH(femmodel->elements,femmodel->parameters);
+
+        /*Solve KU=F*/
         SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
-        _error_("STOP");
         CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
         Reduceloadx(pf, Kfs, ys); delete Kfs;
@@ -30 +37 @@
         delete Kff; delete pf; delete df;
         Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys;
+
+        /*Update solution*/
         InputUpdateFromSolutionx(femmodel,ug); 
+
+        /*Update d and tau accordingly*/
+        _error_("STOP");
+
         delete ug;  
+        delete analysis;
 }