Changeset 17539

Timestamp:

03/25/14 12:09:03 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: working on initialization

Location:

issm/trunk-jpl/src/c/analyses

Files:

• StressbalanceAnalysis.cpp (modified) (2 diffs)
• StressbalanceAnalysis.h (modified) (1 diff)

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){
@@ -820 +821 @@
         }
 }/*}}}*/
+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*/

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*/