Changeset 17109


Ignore:
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:
8 edited

Legend:

Unmodified
Added
Removed

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

    r17108 r17109  
    55#include "../modules/modules.h"
    66#include "../solutionsequences/solutionsequences.h"
     7
     8#define FSANALYTICAL 1
    79
    810/*Model processing*/
    … …  
    27192721ElementVector* StressbalanceAnalysis::CreatePVectorFS(Element* element){/*{{{*/
    27202722
     2723        #ifdef FSANALYTICAL
     2724        ElementVector* pe=CreatePVectorFSAnalytical(element);
     2725
     2726        #else
    27212727        /*compute all load vectors for this element*/
    27222728        ElementVector* pe1=CreatePVectorFSViscous(element);
    … …  
    27292735        delete pe2;
    27302736        delete pe3;
     2737        #endif
     2738
     2739        return pe;
     2740}/*}}}*/
     2741ElementVector* StressbalanceAnalysis::CreatePVectorFSAnalytical(Element* element){/*{{{*/
     2742
     2743        bool        analytical_solution=0;
     2744        int         i,meshtype,dim;
     2745        IssmDouble  x_coord,y_coord,z_coord;
     2746        IssmDouble  Jdet,forcex,forcey,forcez;
     2747        IssmDouble *xyz_list = NULL;
     2748
     2749        element->FindParam(&meshtype,MeshTypeEnum);
     2750        switch(meshtype){
     2751                case Mesh3DEnum:         dim = 3; break;
     2752                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     2753        }
     2754
     2755        /*Fetch number of nodes and dof for this finite element*/
     2756        int vnumnodes = element->GetNumberOfNodesVelocity();
     2757        int pnumnodes = element->GetNumberOfNodesPressure();
     2758
     2759        /*Prepare coordinate system list*/
     2760        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     2761        if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
     2762        else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
     2763        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     2764
     2765        /*Initialize vectors*/
     2766        ElementVector* pe     = element->NewElementVector(FSvelocityEnum);
     2767        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
     2768
     2769        /*Retrieve all inputs and parameters*/
     2770        element->GetVerticesCoordinates(&xyz_list);
     2771        IssmDouble  rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum);
     2772
     2773        /* Start  looping on the number of gaussian points: */
     2774        Gauss* gauss=element->NewGauss(5);
     2775        for(int ig=gauss->begin();ig<gauss->end();ig++){
     2776                gauss->GaussPoint(ig);
     2777
     2778                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     2779                element->NodalFunctionsVelocity(vbasis,gauss);
     2780
     2781                forcex=1.0; forcey=1.0; forcez=1.0;
     2782//              forcex->fx(x_coord,y_coord,z_coord);
     2783//              forcex->fy(x_coord,y_coord,z_coord);
     2784//              forcex->fz(x_coord,y_coord,z_coord);
     2785
     2786                for(i=0;i<vnumnodes;i++){
     2787                        pe->values[i*dim+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i];
     2788                        pe->values[i*dim+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i];
     2789                        pe->values[i*dim+2] += +rho_ice*forcez *Jdet*gauss->weight*vbasis[i];
     2790                }
     2791        }
     2792
     2793        /*Transform coordinate system*/
     2794        element->TransformLoadVectorCoord(pe,cs_list);
     2795
     2796        /*Clean up and return*/
     2797        delete gauss;
     2798        xDelete<int>(cs_list);
     2799        xDelete<IssmDouble>(vbasis);
     2800        xDelete<IssmDouble>(xyz_list);
    27312801        return pe;
    27322802}/*}}}*/

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

    r17075 r17109  
    6969                ElementMatrix* CreateKMatrixFSFriction(Element* element);
    7070                ElementVector* CreatePVectorFS(Element* element);
     71                ElementVector* CreatePVectorFSAnalytical(Element* element);
    7172                ElementVector* CreatePVectorFSViscous(Element* element);
    7273                ElementVector* CreatePVectorFSShelf(Element* element);

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

    r17100 r17109  
    166166                virtual void   GetVerticesCoordinatesBase(IssmDouble** xyz_list)=0;
    167167                virtual void   GetVerticesCoordinatesTop(IssmDouble** xyz_list)=0;
     168                virtual IssmDouble GetXcoord(Gauss* gauss)=0;
     169                virtual IssmDouble GetYcoord(Gauss* gauss)=0;
    168170                virtual IssmDouble GetZcoord(Gauss* gauss)=0;
    169                 virtual IssmDouble GetYcoord(Gauss* gauss)=0;
    170171                virtual void   GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0;
    171172                virtual int    GetElementType(void)=0;

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

    r17065 r17109  
    12071207        /*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
    12081208        input->GetVectorFromInputs(vector,&vertexpidlist[0]);
     1209}
     1210/*}}}*/
     1211/*FUNCTION Penta::GetXcoord {{{*/
     1212IssmDouble Penta::GetXcoord(Gauss* gauss){
     1213
     1214        int    i;
     1215        IssmDouble x;
     1216        IssmDouble xyz_list[NUMVERTICES][3];
     1217        IssmDouble x_list[NUMVERTICES];
     1218
     1219        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     1220        for(i=0;i<NUMVERTICES;i++) x_list[i]=xyz_list[i][0];
     1221        PentaRef::GetInputValue(&x,x_list,gauss);
     1222
     1223        return x;
     1224}
     1225/*}}}*/
     1226/*FUNCTION Penta::GetYcoord {{{*/
     1227IssmDouble Penta::GetYcoord(Gauss* gauss){
     1228
     1229        int    i;
     1230        IssmDouble y;
     1231        IssmDouble xyz_list[NUMVERTICES][3];
     1232        IssmDouble y_list[NUMVERTICES];
     1233
     1234        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     1235        for(i=0;i<NUMVERTICES;i++) y_list[i]=xyz_list[i][1];
     1236        PentaRef::GetInputValue(&y,y_list,gauss);
     1237
     1238        return y;
    12091239}
    12101240/*}}}*/

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

    r17047 r17109  
    9595                int    GetNumberOfVertices(void);
    9696                void   GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type);
     97                IssmDouble GetXcoord(Gauss* gauss);
     98                IssmDouble GetYcoord(Gauss* gauss);
    9799                IssmDouble GetZcoord(Gauss* gauss);
    98                 IssmDouble GetYcoord(Gauss* gauss){_error_("Not implemented");};
    99100                void   GetVectorFromInputs(Vector<IssmDouble>* vector,int name_enum);
    100101                void   GetVerticesCoordinates(IssmDouble** pxyz_list);

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

    r17047 r17109  
    139139                void        GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype){_error_("not implemented yet");};
    140140                Node*       GetNode(int node_number){_error_("Not implemented");};
     141                IssmDouble  GetXcoord(Gauss* gauss){_error_("Not implemented");};
    141142                IssmDouble  GetYcoord(Gauss* gauss){_error_("Not implemented");};
    142143                IssmDouble  GetZcoord(Gauss* gauss){_error_("not implemented yet");};

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

    r17047 r17109  
    212212
    213213                void             GetVertexPidList(int* doflist);
     214                IssmDouble     GetXcoord(Gauss* gauss){_error_("not implemented");};
    214215                IssmDouble     GetYcoord(Gauss* gauss);
    215216                IssmDouble     GetZcoord(Gauss* gauss){_error_("not implemented");};

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

    r15306 r17109  
    4747        MatSetFromOptions(this->matrix);
    4848        MatMPIAIJSetPreallocation(this->matrix,0,d_nnz,0,o_nnz);
    49         //MatSetOption(this->matrix,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
     49        MatSetOption(this->matrix,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
    5050
    5151}
Note: See TracChangeset for help on using the changeset viewer.