Changeset 17852

Timestamp:

04/25/14 15:54:36 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: preparing thermal with tets

Location:

issm/trunk-jpl/src/c

Files:

• classes/Elements/Tetra.cpp (modified) (3 diffs)
• classes/Elements/Tetra.h (modified) (4 diffs)
• modules/ModelProcessorx/CreateSingleNodeToElementConnectivity.cpp (modified) (1 diff)

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

@@ -103 +103 @@
 }
 /*}}}*/
+/*FUNCTION Tetra::ElementSizes (THIS ONE){{{*/
+void Tetra::ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){
+
+        IssmDouble xyz_list[NUMVERTICES][3];
+        IssmDouble xmin,ymin,zmin;
+        IssmDouble xmax,ymax,zmax;
+
+        /*Get xyz list: */
+        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+        xmin=xyz_list[0][0]; xmax=xyz_list[0][0];
+        ymin=xyz_list[0][1]; ymax=xyz_list[0][1];
+        zmin=xyz_list[0][2]; zmax=xyz_list[0][2];
+
+        for(int i=1;i<NUMVERTICES;i++){
+                if(xyz_list[i][0]<xmin) xmin=xyz_list[i][0];
+                if(xyz_list[i][0]>xmax) xmax=xyz_list[i][0];
+                if(xyz_list[i][1]<ymin) ymin=xyz_list[i][1];
+                if(xyz_list[i][1]>ymax) ymax=xyz_list[i][1];
+                if(xyz_list[i][2]<zmin) zmin=xyz_list[i][2];
+                if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2];
+        }
+
+        *hx=xmax-xmin;
+        *hy=ymax-ymin;
+        *hz=zmax-zmin;
+}
+/*}}}*/
 /*FUNCTION Tetra::FaceOnBaseIndices{{{*/
 void Tetra::FaceOnBaseIndices(int* pindex1,int* pindex2,int* pindex3){
@@ -171 +198 @@
         /*return TetraRef field*/
         return this->element_type;
+}
+/*}}}*/
+/*FUNCTION Tetra::GetInputValue(IssmDouble* pvalue,Node* node,int enumtype) (THIS ONE??){{{*/
+void Tetra::GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){
+
+        Input* input=inputs->GetInput(enumtype);
+        if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
+
+        GaussTetra* gauss=new GaussTetra();
+        gauss->GaussVertex(this->GetNodeIndex(node));
+
+        input->GetInputValue(pvalue,gauss);
+        delete gauss;
+}
+/*}}}*/
+/*FUNCTION Tetra::GetNodeIndex (THIS ONE){{{*/
+int Tetra::GetNodeIndex(Node* node){
+
+        _assert_(nodes);
+        int numnodes = this->NumberofNodes();
+
+        for(int i=0;i<numnodes;i++){
+                if(node==nodes[i]) return i;
+        }
+        _error_("Node provided not found among element nodes");
+
 }
 /*}}}*/
@@ -917 +970 @@
 }
 /*}}}*/
+/*FUNCTION Tetra::ViscousHeating (THIS ONE){{{*/
+void Tetra::ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){
+
+        /*Intermediaries*/
+        IssmDouble phi;
+        IssmDouble viscosity;
+        IssmDouble epsilon[6];
+
+        _assert_(gauss->Enum()==GaussTetraEnum);
+        this->StrainRateFS(&epsilon[0],xyz_list,(GaussTetra*)gauss,vx_input,vy_input,vz_input);
+        this->ViscosityFS(&viscosity,3,xyz_list,(GaussTetra*)gauss,vx_input,vy_input,vz_input);
+        GetPhi(&phi,&epsilon[0],viscosity);
+
+        /*Assign output pointer*/
+        *pphi = phi;
+}
+/*}}}*/
 /*FUNCTION Tetra::PressureInterpolation{{{*/
 int Tetra::PressureInterpolation(void){

issm/trunk-jpl/src/c/classes/Elements/Tetra.h

@@ -62 +62 @@
                 void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
                 void        Delta18oParameterization(void){_error_("not implemented yet");};
-                void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");};
+                void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
                 void        FaceOnFrontIndices(int* pindex1,int* pindex2,int* pindex3);
                 void        FaceOnBaseIndices(int* pindex1,int* pindex2,int* pindex3);
@@ -72 +72 @@
                 Element*    GetSurfaceElement(void){_error_("not implemented yet");};
                 Element*    GetBasalElement(void){_error_("not implemented yet");};
-                int         GetNodeIndex(Node* node){_error_("not implemented yet");};
+                int         GetNodeIndex(Node* node);
                 int         GetNumberOfNodes(void);
                 int         GetNumberOfVertices(void);
@@ -112 +112 @@
                 void        GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating){_error_("not implemented yet");};
                 IssmDouble  GetGroundedPortion(IssmDouble* xyz_list){_error_("not implemented yet");};
-                void        GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){_error_("not implemented yet");};
+                void        GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
                 Node*       GetNode(int node_number){_error_("Not implemented");};
                 IssmDouble  GetXcoord(Gauss* gauss);
@@ -128 +128 @@
                 int         VertexConnectivity(int vertexindex){_error_("not implemented yet");};
                 void        VerticalSegmentIndices(int** pindices,int* pnumseg){_error_("not implemented yet");};
-                void        ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented yet");};
+                void        ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
                 bool        IsZeroLevelset(int levelset_enum){_error_("not implemented");};
                 bool               IsIcefront(void);

issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateSingleNodeToElementConnectivity.cpp

@@ -29 +29 @@
         connectivity=xNewZeroInit<int>(iomodel->numberofvertices);
 
-        /*Get element width (3 or 6)*/
-        switch(iomodel->domaintype){
-                case Domain2DhorizontalEnum: elementswidth=3; break;
-                case Domain2DverticalEnum  : elementswidth=3; break;
-                case Domain3DEnum          : elementswidth=6; break;
+        /*Get element width*/
+        switch(iomodel->meshelementtype){
+                case TriaEnum:  elementswidth=3; break;
+                case PentaEnum: elementswidth=6; break;
+                case TetraEnum: elementswidth=4; break;
                 default:  _error_("mesh type "<<EnumToStringx(iomodel->domaintype)<<" not supported yet");
         }