Changeset 17288

Timestamp:

02/15/14 20:33:24 (11 years ago)

Author:

seroussi

Message:

BUG: fixing analytical 2d

Location:

issm/trunk-jpl/src/c

Files:

• analyses/StressbalanceAnalysis.cpp (modified) (2 diffs)
• classes/Elements/Tria.cpp (modified) (1 diff)
• classes/Elements/Tria.h (modified) (1 diff)
• shared/FSanalyticals/FSanalytical3.cpp (modified) (1 diff)

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

@@ -6 +6 @@
 #include "../solutionsequences/solutionsequences.h"
 
-#define FSANALYTICAL 3
+//#define FSANALYTICAL 3
 
 /*Model processing*/
@@ -409 +409 @@
                                         if(nodeonbed[i]>0. && groundedice_ls[i]>0. && nodeonFS[i]>0.){
                                                 if(vertices_type[i] == FSApproximationEnum){
-                                                //      for(j=0;j<Nz;j++) spcvz[i*Nz+j] = 0.;
+                                                        for(j=0;j<Nz;j++) spcvz[i*Nz+j] = 0.;
                                                 }
                                                 else{

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

@@ -1191 +1191 @@
         /*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
         input->GetVectorFromInputs(vector,&vertexpidlist[0]);
+}
+/*}}}*/
+/*FUNCTION Tria::GetXcoord {{{*/
+IssmDouble Tria::GetXcoord(Gauss* gauss){
+
+        IssmDouble x;
+        IssmDouble xyz_list[NUMVERTICES][3];
+        IssmDouble x_list[NUMVERTICES];
+
+        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+        for(int i=0;i<NUMVERTICES;i++) x_list[i]=xyz_list[i][0];
+        TriaRef::GetInputValue(&x,x_list,gauss,P1Enum);
+
+        return x;
 }
 /*}}}*/

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

@@ -209 +209 @@
 
                 void             GetVertexPidList(int* doflist);
-                IssmDouble     GetXcoord(Gauss* gauss){_error_("not implemented");};
+                IssmDouble     GetXcoord(Gauss* gauss);
                 IssmDouble     GetYcoord(Gauss* gauss);
                 IssmDouble     GetZcoord(Gauss* gauss){_error_("not implemented");};

issm/trunk-jpl/src/c/shared/FSanalyticals/FSanalytical3.cpp

@@ -8 +8 @@
 #include "../Numerics/constants.h"
 
-IssmDouble fx2(IssmDouble x,IssmDouble y{   
-        IssmDouble p = 2.0;  
-        IssmDouble q = 2.0;   
+IssmDouble fx3(IssmDouble x,IssmDouble y){   
         IssmDouble functionx;
+        IssmDouble z=y;
 
-        functionx =  -2*(y-1)*(y-1)*(x-1)*(2*x-1)*y-4*(y-1)*(y-1)*(x-1)*x*y-2*(y-1)*(y-1)*(2*x-1)*x*y +4*(y-1)*(2*y-1)*(x-1)*(x-1)*y+16*(y-1)*(2*y-1)*(x-1)*x*y+4*(y-1)*(2*y-1)*x*x*y +4*(y-1)*(x-1)*(x-1)*x*x-2*(y-1)*(x-1)*(2*x-1)*y*y-4*(y-1)*(x-1)*x*y*y -2*(y-1)*(2*x-1)*x*y*y+2*(2*y-1)*(x-1)*(x-1)*x*x+4*(x-1)*(x-1)*x*x*y+y-1/2; 
+        functionx =  4*pow(x, 2)*z*pow(x - 1, 2) + 4*pow(x, 2)*z*(z - 1)*(2*z - 1) + 4*pow(x, 2)*pow(x - 1, 2)*(z - 1) + 2*pow(x, 2)*pow(x - 1, 2)*(2*z - 1) + 16*x*z*(x - 1)*(z - 1)*(2*z - 1) - 4*pow(z, 3)*(6*pow(x, 2) - 6*x + 1) + 6*pow(z, 2)*(6*pow(x, 2) - 6*x + 1) + 4*z*pow(x - 1, 2)*(z - 1)*(2*z - 1) - 2*z*(6*pow(x, 2) - 6*x + 1) + z - 1.0L/2.0L;
         return functionx;
 }
-IssmDouble fy2(IssmDouble x,IssmDouble y){   
-        IssmDouble p = 2.0;  
-        IssmDouble q = 2.0;   
+IssmDouble fy3(IssmDouble x,IssmDouble y){   
         IssmDouble functiony;
+        IssmDouble z=y;
 
-        functiony = -4*(y-1)*(y-1)*(x-1)*(2*x-1)*x-4*(y-1)*(y-1)*(x-1)*y*y-2*(y-1)*(y-1)*(2*x-1)*y*y -4*(y-1)*(y-1)*x*y*y+2*(y-1)*(2*y-1)*(x-1)*(x-1)*x+2*(y-1)*(2*y-1)*(x-1)*x*x +4*(y-1)*(x-1)*(x-1)*x*y-16*(y-1)*(x-1)*(2*x-1)*x*y+4*(y-1)*(x-1)*x*x*y +2*(2*y-1)*(x-1)*(x-1)*x*y+2*(2*y-1)*(x-1)*x*x*y -4*(x-1)*(2*x-1)*x*y*y+x-1/2;
+        functiony = -8*pow(x, 3) + 4*pow(x, 2)*z*(x - 1)*(z - 1) + 2*pow(x, 2)*z*(x - 1)*(2*z - 1) + 2*pow(x, 2)*(x - 1)*(z - 1)*(2*z - 1) + 12*pow(x, 2) + 4*x*z*pow(x - 1, 2)*(z - 1) + 2*x*z*pow(x - 1, 2)*(2*z - 1) + 2*x*pow(x - 1, 2)*(z - 1)*(2*z - 1) - 3*x - 6*pow(z, 4)*(2*x - 1) + 12*pow(z, 3)*(2*x - 1) - 6*pow(z, 2)*(2*x - 1) - 24*pow(z, 2)*(2*pow(x, 3) - 3*pow(x, 2) + x) + 24*z*(2*pow(x, 3) - 3*pow(x, 2) + x) - 1.0L/2.0L;
 
         return functiony;