Changeset 17109

Timestamp:

01/14/14 13:15:26 (11 years ago)

Author:

seroussi

Message:

NEW: starting to add analytical functions for FS

Location:

issm/trunk-jpl/src/c

Files:

• analyses/StressbalanceAnalysis.cpp (modified) (3 diffs)
• analyses/StressbalanceAnalysis.h (modified) (1 diff)
• classes/Elements/Element.h (modified) (1 diff)
• classes/Elements/Penta.cpp (modified) (1 diff)
• classes/Elements/Penta.h (modified) (1 diff)
• classes/Elements/Seg.h (modified) (1 diff)
• classes/Elements/Tria.h (modified) (1 diff)
• toolkits/petsc/objects/PetscMat.cpp (modified) (1 diff)

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

@@ -5 +5 @@
 #include "../modules/modules.h"
 #include "../solutionsequences/solutionsequences.h"
+
+#define FSANALYTICAL 1
 
 /*Model processing*/
@@ -2719 +2721 @@
 ElementVector* StressbalanceAnalysis::CreatePVectorFS(Element* element){/*{{{*/
 
+        #ifdef FSANALYTICAL
+        ElementVector* pe=CreatePVectorFSAnalytical(element);
+
+        #else
         /*compute all load vectors for this element*/
         ElementVector* pe1=CreatePVectorFSViscous(element);
@@ -2729 +2735 @@
         delete pe2;
         delete pe3;
+        #endif
+
+        return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorFSAnalytical(Element* element){/*{{{*/
+
+        bool        analytical_solution=0;
+        int         i,meshtype,dim;
+        IssmDouble  x_coord,y_coord,z_coord;
+        IssmDouble  Jdet,forcex,forcey,forcez;
+        IssmDouble *xyz_list = NULL;
+
+        element->FindParam(&meshtype,MeshTypeEnum);
+        switch(meshtype){
+                case Mesh3DEnum:         dim = 3; break;
+                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+        }
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int vnumnodes = element->GetNumberOfNodesVelocity();
+        int pnumnodes = element->GetNumberOfNodesPressure();
+
+        /*Prepare coordinate system list*/
+        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
+        if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
+        else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
+        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
+
+        /*Initialize vectors*/
+        ElementVector* pe     = element->NewElementVector(FSvelocityEnum);
+        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinates(&xyz_list);
+        IssmDouble  rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=element->NewGauss(5);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                element->NodalFunctionsVelocity(vbasis,gauss);
+
+                forcex=1.0; forcey=1.0; forcez=1.0;
+//              forcex->fx(x_coord,y_coord,z_coord);
+//              forcex->fy(x_coord,y_coord,z_coord);
+//              forcex->fz(x_coord,y_coord,z_coord);
+
+                for(i=0;i<vnumnodes;i++){
+                        pe->values[i*dim+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i];
+                        pe->values[i*dim+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i];
+                        pe->values[i*dim+2] += +rho_ice*forcez *Jdet*gauss->weight*vbasis[i];
+                }
+        }
+
+        /*Transform coordinate system*/
+        element->TransformLoadVectorCoord(pe,cs_list);
+
+        /*Clean up and return*/
+        delete gauss;
+        xDelete<int>(cs_list);
+        xDelete<IssmDouble>(vbasis);
+        xDelete<IssmDouble>(xyz_list);
         return pe;
 }/*}}}*/

issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h ¶

@@ -69 +69 @@
                 ElementMatrix* CreateKMatrixFSFriction(Element* element);
                 ElementVector* CreatePVectorFS(Element* element);
+                ElementVector* CreatePVectorFSAnalytical(Element* element);
                 ElementVector* CreatePVectorFSViscous(Element* element);
                 ElementVector* CreatePVectorFSShelf(Element* element);

issm/trunk-jpl/src/c/classes/Elements/Element.h ¶

@@ -166 +166 @@
                 virtual void   GetVerticesCoordinatesBase(IssmDouble** xyz_list)=0;
                 virtual void   GetVerticesCoordinatesTop(IssmDouble** xyz_list)=0;
+                virtual IssmDouble GetXcoord(Gauss* gauss)=0;
+                virtual IssmDouble GetYcoord(Gauss* gauss)=0;
                 virtual IssmDouble GetZcoord(Gauss* gauss)=0;
-                virtual IssmDouble GetYcoord(Gauss* gauss)=0;
                 virtual void   GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0;
                 virtual int    GetElementType(void)=0;

issm/trunk-jpl/src/c/classes/Elements/Penta.cpp ¶

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

issm/trunk-jpl/src/c/classes/Elements/Penta.h ¶

@@ -95 +95 @@
                 int    GetNumberOfVertices(void);
                 void   GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type);
+                IssmDouble GetXcoord(Gauss* gauss);
+                IssmDouble GetYcoord(Gauss* gauss);
                 IssmDouble GetZcoord(Gauss* gauss);
-                IssmDouble GetYcoord(Gauss* gauss){_error_("Not implemented");};
                 void   GetVectorFromInputs(Vector<IssmDouble>* vector,int name_enum);
                 void   GetVerticesCoordinates(IssmDouble** pxyz_list);

issm/trunk-jpl/src/c/classes/Elements/Seg.h ¶

@@ -139 +139 @@
                 void        GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype){_error_("not implemented yet");};
                 Node*       GetNode(int node_number){_error_("Not implemented");};
+                IssmDouble  GetXcoord(Gauss* gauss){_error_("Not implemented");};
                 IssmDouble  GetYcoord(Gauss* gauss){_error_("Not implemented");};
                 IssmDouble  GetZcoord(Gauss* gauss){_error_("not implemented yet");};

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

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

issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscMat.cpp ¶

@@ -47 +47 @@
         MatSetFromOptions(this->matrix);
         MatMPIAIJSetPreallocation(this->matrix,0,d_nnz,0,o_nnz);
-        //MatSetOption(this->matrix,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
+        MatSetOption(this->matrix,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
 
 }