Changeset 350

Timestamp:

05/12/09 10:43:05 (16 years ago)

Author:

seroussi

Message:

added pressure comptute in sing and beam elements

Location:

issm/trunk/src/c/objects

Files:

• Beam.cpp (modified) (1 diff)
• Penta.cpp (modified) (1 diff)
• Sing.cpp (modified) (1 diff)

issm/trunk/src/c/objects/Beam.cpp

@@ -565 +565 @@
 
 #undef __FUNCT__ 
+#define __FUNCT__ "Beam::GetParameterValue"
+void Beam::GetParameterValue(double* pvalue, double* value_list,double gauss_coord){
+
+        double l1l2[2];
+        
+        GetNodalFunctions(&l1l2[0],gauss_coord);
+
+        *pvalue=l1l2[0]*value_list[0]+l1l2[1]*value_list[1];
+}
+
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Beam::GetNodalFunctions"
+void Beam::GetNodalFunctions(double* l1l2, double gauss_coord){
+        
+        /*This routine returns the values of the nodal functions  at the gaussian point.*/
+
+        /*First nodal function: */
+        l1l2[0]=.5*gauss_coord+.5;
+
+        /*Second nodal function: */
+        l1l2[1]=-.5*gauss_coord+.5;
+}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Beam::GetJacobianDeterminant"
+void Beam::GetJacobianDeterminant(double* pJdet,double* z_list, double gauss_coord){
+
+
+        double Jdet;
+
+        Jdet=1.0/2.0*(z_list[1]-z_list[0]);
+
+        if(Jdet<0){
+                throw ErrorException(__FUNCT__," negative jacobian determinant!");
+        }
+        
+        *pJdet=Jdet;
+}
+
+#undef __FUNCT__ 
 #define __FUNCT__ "Beam::ComputePressure"
 void  Beam::ComputePressure(Vec p_g){
-        throw ErrorException(__FUNCT__," not supported yet!");
-}
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Beam::GetParameterValue"
-void Beam::GetParameterValue(double* pvalue, double* value_list,double gauss_coord){
-
-        double l1l2[2];
-        
-        GetNodalFunctions(&l1l2[0],gauss_coord);
-
-        *pvalue=l1l2[0]*value_list[0]+l1l2[1]*value_list[1];
-}
-
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Beam::GetNodalFunctions"
-void Beam::GetNodalFunctions(double* l1l2, double gauss_coord){
-        
-        /*This routine returns the values of the nodal functions  at the gaussian point.*/
-
-        /*First nodal function: */
-        l1l2[0]=.5*gauss_coord+.5;
-
-        /*Second nodal function: */
-        l1l2[1]=-.5*gauss_coord+.5;
-}
-
-#undef __FUNCT__ 
-#define __FUNCT__ "Beam::GetJacobianDeterminant"
-void Beam::GetJacobianDeterminant(double* pJdet,double* z_list, double gauss_coord){
-
-
-        double Jdet;
-
-        Jdet=1.0/2.0*(z_list[1]-z_list[0]);
-
-        if(Jdet<0){
-                throw ErrorException(__FUNCT__," negative jacobian determinant!");
-        }
-        
-        *pJdet=Jdet;
-}
-
+
+        int i;
+        const int numgrids=2;
+        int doflist[numgrids];
+        double pressure[numgrids];
+        double rho_ice,g;
+        double xyz_list[numgrids][3];
+
+        /*Get node data: */
+        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+        
+        /*Get dof list on which we will plug the pressure values: */
+        GetDofList1(&doflist[0]);
+
+        /*pressure is lithostatic: */
+        rho_ice=matpar->GetRhoIce();
+        g=matpar->GetG();
+        for(i=0;i<numgrids;i++){
+                pressure[i]=rho_ice*g*(s[i]-xyz_list[i][2]);
+        }
+        
+        /*plug local pressure values into global pressure vector: */
+        VecSetValues(p_g,numgrids,doflist,(const double*)pressure,INSERT_VALUES);
+
+}
+

issm/trunk/src/c/objects/Penta.cpp

@@ -1770 +1770 @@
         double pressure[numgrids];
         double rho_ice,g;
-        double       xyz_list[numgrids][3];
+        double xyz_list[numgrids][3];
                 
         /*Get node data: */

issm/trunk/src/c/objects/Sing.cpp

@@ -460 +460 @@
 #define __FUNCT__ "Sing::ComputePressure"
 void  Sing::ComputePressure(Vec p_g){
-        throw ErrorException(__FUNCT__," not supported yet!");
-}
+
+        int i;
+        const int numgrids=1;
+        int doflist[numgrids];
+        double pressure[numgrids];
+        double rho_ice,g;
+
+        /*Get dof list on which we will plug the pressure values: */
+        GetDofList1(&doflist[0]);
+
+        /*pressure is lithostatic: */
+        rho_ice=matpar->GetRhoIce();
+        g=matpar->GetG();
+        pressure[0]=rho_ice*g*h;
+        
+        /*plug local pressure values into global pressure vector: */
+        VecSetValues(p_g,numgrids,doflist,(const double*)pressure,INSERT_VALUES);
+
+}