Changeset 3800

Timestamp:

05/18/10 10:47:41 (15 years ago)

Author:

Mathieu Morlighem

Message:

minor

Location:

issm/trunk/src/c

Files:

• DataSet/Inputs.cpp (modified) (1 diff)
• objects/Elements/Tria.cpp (modified) (1 diff)
• objects/Inputs/TriaVertexInput.cpp (modified) (3 diffs)
• objects/Inputs/TriaVertexInput.h (modified) (1 diff)
• parallel/prognostic_core.cpp (modified) (1 diff)

issm/trunk/src/c/DataSet/Inputs.cpp

@@ -300 +300 @@
         }
 
-        if (!xinput | !yinput){
+        if (!xinput || !yinput){
                 /*we could not find one input with the correct enum type. No defaults values were provided, 
                  * error out: */

issm/trunk/src/c/objects/Elements/Tria.cpp

@@ -1211 +1211 @@
                 /*Compute thickness at gaussian point: */
                 inputs->GetParameterValue(&thickness, gauss_l1l2l3,ThicknessEnum);
+                printf("thickness = %g\n",thickness);
 
                 /*Get strain rate from velocity: */

issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp

@@ -149 +149 @@
 /*}}}*/
 /*FUNCTION TriaVertexInput::GetParameterValue(double* pvalue,double* gauss){{{1*/
-void TriaVertexInput::GetParameterValue(double* pvalue,double* gauss){ISSMERROR(" not supported yet!");}
+void TriaVertexInput::GetParameterValue(double* pvalue,double* gauss){
+        /*P1 interpolation on Gauss point*/
+
+        /*intermediary*/
+        double l1l2l3[3];
+
+        /*nodal functions: */
+        GetNodalFunctions(l1l2l3,gauss);
+
+        /*Assign output pointers:*/
+        *pvalue=l1l2l3[0]*values[0]+l1l2l3[1]*values[1]+l1l2l3[2]*values[2];
+
+}
 /*}}}*/
 /*FUNCTION TriaVertexInput::GetParameterValue(double* pvalue,double* gauss,double defaultvalue){{{1*/
@@ -161 +173 @@
 /*}}}*/
 /*FUNCTION TriaVertexInput::GetStrainRate(double* epsilon,Input* yinput, double* xyz_list, double* gauss){{{1*/
-void TriaVertexInput::GetStrainRate(double* epsilon,Input* yinput,double* xyz_list, double* gauss){ISSMERROR(" not supported yet!");}
+void TriaVertexInput::GetStrainRate(double* epsilon,Input* yinput,double* xyz_list, double* gauss){
+        ISSMERROR("not implemented yet");
+}
 /*}}}*/
 /*FUNCTION TriaVertexInput::GetStrainRateStokes(double* epsilon,Input* yinput, Input* zinput, double* xyz_list, double* gauss){{{1*/
@@ -171 +185 @@
 }
 /*}}}*/
+
+/*Intermediary*/
+/*FUNCTION TriaVertexInput::GetNodalFunctions {{{1*/
+void TriaVertexInput::GetNodalFunctions(double* l1l2l3, double* gauss){
+        
+        /*This routine returns the values of the nodal functions  at the gaussian point.*/
+
+        /*First nodal function: */
+        l1l2l3[0]=gauss[0];
+
+        /*Second nodal function: */
+        l1l2l3[1]=gauss[1];
+
+        /*Third nodal function: */
+        l1l2l3[2]=gauss[2];
+
+}
+/*}}}*/
+/*FUNCTION TriaVertexInput::GetB {{{1*/
+void TriaVertexInput::GetB(double* B, double* xyz_list, double* gauss){
+        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
+         * For grid i, Bi can be expressed in the actual coordinate system
+         * by: 
+         *       Bi=[ dh/dx           0    ]
+         *          [   0           dh/dy  ]
+         *          [ 1/2*dh/dy  1/2*dh/dx ]
+         * where h is the interpolation function for grid i.
+         *
+         * We assume B has been allocated already, of size: 3x(NDOF2*numgrids)
+         */
+
+        int i;
+        const int NDOF2=2;
+        const int numgrids=3;
+
+        double dh1dh3[NDOF2][numgrids];
+
+
+        /*Get dh1dh2dh3 in actual coordinate system: */
+        GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);
+
+        /*Build B: */
+        for (i=0;i<numgrids;i++){
+                *(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh3[0][i]; //B[0][NDOF2*i]=dh1dh3[0][i];
+                *(B+NDOF2*numgrids*0+NDOF2*i+1)=0;
+                *(B+NDOF2*numgrids*1+NDOF2*i)=0;
+                *(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh3[1][i];
+                *(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh3[1][i]; 
+                *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh3[0][i]; 
+        }
+}
+/*}}}*/
+/*FUNCTION TriaVertexInput::GetNodalFunctionsDerivatives {{{1*/
+void TriaVertexInput::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss){
+
+        /*This routine returns the values of the nodal functions derivatives  (with respect to the 
+         * actual coordinate system: */
+
+        int i;
+        const int NDOF2=2;
+        const int numgrids=3;
+
+        double dh1dh3_ref[NDOF2][numgrids];
+        double Jinv[NDOF2][NDOF2];
+
+
+        /*Get derivative values with respect to parametric coordinate system: */
+        GetNodalFunctionsDerivativesReference(&dh1dh3_ref[0][0], gauss); 
+
+        /*Get Jacobian invert: */
+        GetJacobianInvert(&Jinv[0][0], xyz_list, gauss);
+
+        /*Build dh1dh3: 
+         *
+         * [dhi/dx]= Jinv*[dhi/dr]
+         * [dhi/dy]       [dhi/ds]
+         */
+
+        for (i=0;i<numgrids;i++){
+                *(dh1dh3+numgrids*0+i)=Jinv[0][0]*dh1dh3_ref[0][i]+Jinv[0][1]*dh1dh3_ref[1][i];
+                *(dh1dh3+numgrids*1+i)=Jinv[1][0]*dh1dh3_ref[0][i]+Jinv[1][1]*dh1dh3_ref[1][i];
+        }
+
+}
+/*}}}*/
+/*FUNCTION TriaVertexInput::GetNodalFunctionsDerivativesReference {{{1*/
+void TriaVertexInput::GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss_l1l2l3){
+
+        /*This routine returns the values of the nodal functions derivatives  (with respect to the 
+         * natural coordinate system) at the gaussian point. */
+
+        const int NDOF2=2;
+        const int numgrids=3;
+
+        /*First nodal function: */
+        *(dl1dl3+numgrids*0+0)=-0.5; 
+        *(dl1dl3+numgrids*1+0)=-1.0/(2.0*SQRT3);
+
+        /*Second nodal function: */
+        *(dl1dl3+numgrids*0+1)=0.5;
+        *(dl1dl3+numgrids*1+1)=-1.0/(2.0*SQRT3);
+
+        /*Third nodal function: */
+        *(dl1dl3+numgrids*0+2)=0;
+        *(dl1dl3+numgrids*1+2)=1.0/SQRT3;
+
+}
+/*}}}*/
+/*FUNCTION TriaVertexInput::GetJacobian {{{1*/
+void TriaVertexInput::GetJacobian(double* J, double* xyz_list,double* gauss){
+
+        /*The Jacobian is constant over the element, discard the gaussian points. 
+         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
+
+        const int NDOF2=2;
+        const int numgrids=3;
+        double x1,y1,x2,y2,x3,y3;
+
+        x1=*(xyz_list+numgrids*0+0);
+        y1=*(xyz_list+numgrids*0+1);
+        x2=*(xyz_list+numgrids*1+0);
+        y2=*(xyz_list+numgrids*1+1);
+        x3=*(xyz_list+numgrids*2+0);
+        y3=*(xyz_list+numgrids*2+1);
+
+
+        *(J+NDOF2*0+0)=0.5*(x2-x1);
+        *(J+NDOF2*1+0)=SQRT3/6.0*(2*x3-x1-x2);
+        *(J+NDOF2*0+1)=0.5*(y2-y1);
+        *(J+NDOF2*1+1)=SQRT3/6.0*(2*y3-y1-y2);
+}
+/*}}}*/
+/*FUNCTION TriaVertexInput::GetJacobianInvert {{{1*/
+void TriaVertexInput::GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss){
+
+        double Jdet;
+        const int NDOF2=2;
+        const int numgrids=3;
+
+        /*Call Jacobian routine to get the jacobian:*/
+        GetJacobian(Jinv, xyz_list, gauss);
+
+        /*Invert Jacobian matrix: */
+        MatrixInverse(Jinv,NDOF2,NDOF2,NULL,0,&Jdet);
+
+}
+/*}}}*/

issm/trunk/src/c/objects/Inputs/TriaVertexInput.h

@@ -63 +63 @@
                 void GetStrainRateStokes(double* epsilon,Input* yinput, Input* zinput, double* xyz_list, double* gauss);
                 void ChangeEnum(int newenumtype);
+
+                void GetNodalFunctions(double* l1l2l3, double* gauss);
+                void GetB(double* B, double* xyz_list, double* gauss);
+                void GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss);
+                void GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss_l1l2l3);
+                void GetJacobian(double* J, double* xyz_list,double* gauss);
+                void GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss);
                 /*}}}*/
 

issm/trunk/src/c/parallel/prognostic_core.cpp

@@ -47 +47 @@
         
         _printf_("extract result from extruded inputs: \n");
-        InputToResult(&result,fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,ThicknessEnum,results->Size()+1,0,1);
+        result=new Result(results->Size()+1,0,1,"h_g",h_g);
+        //InputToResult(&result,fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,ThicknessEnum,results->Size()+1,0,1);
         results->AddObject(result);