Changeset 17359

Timestamp:

02/26/14 20:46:50 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added HO convergence tests

Location:

issm/trunk-jpl/src/c

Files:

• analyses/StressbalanceAnalysis.cpp (modified) (3 diffs)
• shared/FSanalyticals/fsanalyticals.cpp (modified) (3 diffs)
• shared/FSanalyticals/fsanalyticals.h (modified) (1 diff)

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

@@ -7 +7 @@
 #include "../cores/cores.h"
 
-//#define FSANALYTICAL 3
+//#define FSANALYTICAL 101
 
 /*Model processing*/
@@ -2273 +2273 @@
         return Ke;
 }/*}}}*/
+#ifdef FSANALYTICAL
+ElementVector* StressbalanceAnalysis::CreatePVectorHO(Element* element){/*{{{*/
+
+        /*Intermediaries */
+        int         dim,meshtype;
+        IssmDouble  x_coord,y_coord,z_coord;
+        IssmDouble  Jdet,forcex,forcey,forcez;
+        IssmDouble* xyz_list = NULL;
+
+        /*Get problem dimension*/
+        element->FindParam(&meshtype,MeshTypeEnum);
+        switch(meshtype){
+                case Mesh2DverticalEnum: dim = 2; break;
+                case Mesh3DEnum:         dim = 3; break;
+                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+        }
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Initialize Element vector and vectors*/
+        ElementVector* pe=element->NewElementVector(HOApproximationEnum);
+        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinates(&xyz_list);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=element->NewGauss(3);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                element->NodalFunctions(basis, gauss);
+
+                x_coord=element->GetXcoord(gauss);
+                y_coord=element->GetYcoord(gauss);
+                if(dim==3) z_coord=element->GetZcoord(gauss);
+                else z_coord=0.;
+
+                forcex=fx(x_coord,y_coord,z_coord,FSANALYTICAL);
+                forcey=fy(x_coord,y_coord,z_coord,FSANALYTICAL);
+
+                for(int i=0;i<numnodes;i++){
+                        pe->values[i*(dim-1)+0]+=forcex*Jdet*gauss->weight*basis[i];
+                        pe->values[i*(dim-1)+1]+=forcey*Jdet*gauss->weight*basis[i];
+                }
+        }
+
+        /*Transform coordinate system*/
+        if(dim==3) element->TransformLoadVectorCoord(pe,XYEnum);
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(basis);
+        xDelete<IssmDouble>(xyz_list);
+        delete gauss;
+        return pe;
+}/*}}}*/
+#else
 ElementVector* StressbalanceAnalysis::CreatePVectorHO(Element* element){/*{{{*/
 
@@ -2285 +2344 @@
         return pe;
 }/*}}}*/
+#endif
 ElementVector* StressbalanceAnalysis::CreatePVectorHODrivingStress(Element* element){/*{{{*/
 

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

@@ -36 +36 @@
                 case 8: 
                         return fx=1.0;
+
+                case 101: 
+                        return fx=fx101(x,y,z);
                 default:
                         _error_("FS analytical solution"<<testid<<" not implemented yet");
@@ -61 +64 @@
                 case 8: 
                         return fy=1.0;
+
+                case 101: 
+                        return fy=fy101(x,y,z);
                 default:
                         _error_("FS analytical solution"<<testid<<" not implemented yet");
@@ -253 +259 @@
 }
 /*}}}*/
+
+IssmDouble fx101(IssmDouble x,IssmDouble y,IssmDouble z){   /*{{{*/
+
+   return 4*pow(x, 2)*y*z*pow(x - 1, 2)*(z - 1) + 8*pow(x, 2)*y*z*(y - 1)*(2*y - 1)*(z - 1) + 2*pow(x, 2)*y*pow(x - 1, 2)*(y - 1)*(2*y - 1) + 4*pow(x, 2)*z*pow(x - 1, 2)*(y - 1)*(z - 1) + 2*pow(x, 2)*z*pow(x - 1, 2)*(2*y - 1)*(z - 1) - 12*x*pow(y, 2)*z*(x - 1)*(y - 1)*(z - 1) - 6*x*pow(y, 2)*z*(2*x - 1)*(y - 1)*(z - 1) - 12*x*y*z*(x - 1)*pow(y - 1, 2)*(z - 1) + 32*x*y*z*(x - 1)*(y - 1)*(2*y - 1)*(z - 1) - 6*x*y*z*(2*x - 1)*pow(y - 1, 2)*(z - 1) - 6*pow(y, 2)*z*(x - 1)*(2*x - 1)*(y - 1)*(z - 1) + 8*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1) - 6*y*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1);
+}
+/*}}}*/
+IssmDouble fy101(IssmDouble x,IssmDouble y,IssmDouble z){   /*{{{*/
+
+        return 12*pow(x, 2)*y*z*(x - 1)*(y - 1)*(z - 1) + 6*pow(x, 2)*y*z*(x - 1)*(2*y - 1)*(z - 1) + 6*pow(x, 2)*z*(x - 1)*(y - 1)*(2*y - 1)*(z - 1) - 8*x*pow(y, 2)*z*(x - 1)*(2*x - 1)*(z - 1) - 4*x*pow(y, 2)*z*pow(y - 1, 2)*(z - 1) - 2*x*pow(y, 2)*(x - 1)*(2*x - 1)*pow(y - 1, 2) + 12*x*y*z*pow(x - 1, 2)*(y - 1)*(z - 1) + 6*x*y*z*pow(x - 1, 2)*(2*y - 1)*(z - 1) - 32*x*y*z*(x - 1)*(2*x - 1)*(y - 1)*(z - 1) + 6*x*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1) - 8*x*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1) - 4*pow(y, 2)*z*(x - 1)*pow(y - 1, 2)*(z - 1) - 2*pow(y, 2)*z*(2*x - 1)*pow(y - 1, 2)*(z - 1);
+}
+/*}}}*/

issm/trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.h

@@ -30 +30 @@
 IssmDouble fy7(IssmDouble x_coord, IssmDouble y_coord);
 
+IssmDouble fx101(IssmDouble x_coord, IssmDouble y_coord, IssmDouble z_coord);
+IssmDouble fy101(IssmDouble x_coord, IssmDouble y_coord, IssmDouble z_coord);
 #endif //ifndef _SHARED_ANALYTICALS_H_