Changeset 16839


Ignore:
Timestamp:
11/20/13 19:08:46 (11 years ago)
Author:
seroussi
Message:

CHG: moved ice front HO in analysis

Location:
issm/trunk-jpl/src/c
Files:
7 edited

Legend:

Unmodified
Added
Removed

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

    r16838 r16839  
    11361136        element->GetVerticesCoordinates(&xyz_list);
    11371137        element->ZeroLevelsetCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
    1138         element->NormalBase(&normal[0],xyz_list_front);
     1138        element->NormalSection(&normal[0],xyz_list_front);
    11391139
    11401140        /*Start looping on Gaussian points*/
    … …  
    15701570ElementVector* StressbalanceAnalysis::CreatePVectorHOFront(Element* element){/*{{{*/
    15711571
    1572         return NULL;
     1572        /*If no front, return NULL*/
     1573        if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) return NULL;
     1574
     1575        /*Intermediaries*/
     1576        IssmDouble  rho_ice,rho_water,gravity;
     1577        IssmDouble  Jdet,surface,z_g,water_pressure,ice_pressure;
     1578        IssmDouble  surface_under_water,base_under_water,pressure;
     1579        IssmDouble* xyz_list       = NULL;
     1580        IssmDouble* xyz_list_front = NULL;
     1581        IssmDouble  normal[3];
     1582        Gauss*      gauss = NULL;
     1583
     1584        /*Fetch number of nodes and dof for this finite element*/
     1585        int numnodes = element->GetNumberOfNodes();
     1586
     1587        /*Initialize Element vector and other vectors*/
     1588        ElementVector* pe    = element->NewElementVector(HOApproximationEnum);
     1589        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     1590
     1591        /*Retrieve all inputs and parameters*/
     1592        Input* surface_input=element->GetInput(SurfaceEnum); _assert_(surface_input);
     1593        rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum);
     1594        rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
     1595        gravity   = element->GetMaterialParameter(ConstantsGEnum);
     1596        element->GetVerticesCoordinates(&xyz_list);
     1597        element->ZeroLevelsetCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
     1598        element->NormalSection(&normal[0],xyz_list_front);
     1599
     1600        /*Initialize gauss points*/
     1601        IssmDouble zmax=xyz_list[0*3+2]; for(int i=1;i<6;i++) if(xyz_list[i*3+2]>zmax) zmax=xyz_list[i*3+2];
     1602        IssmDouble zmin=xyz_list[0*3+2]; for(int i=1;i<6;i++) if(xyz_list[i*3+2]<zmin) zmin=xyz_list[i*3+2];
     1603        if(zmax>0 && zmin<0) gauss=element->NewGauss(xyz_list,xyz_list_front,3,10);//refined in vertical because of the sea level discontinuity
     1604        else                 gauss=element->NewGauss(xyz_list,xyz_list_front,3,3);
     1605
     1606        /* Start  looping on the number of gaussian points: */
     1607        for(int ig=gauss->begin();ig<gauss->end();ig++){
     1608
     1609                gauss->GaussPoint(ig);
     1610                surface_input->GetInputValue(&surface,gauss);
     1611                z_g=element->GetZcoord(gauss);
     1612                element->NodalFunctions(basis,gauss);
     1613                element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
     1614
     1615                water_pressure = rho_water*gravity*min(0.,z_g);//0 if the gaussian point is above water level
     1616                ice_pressure   = rho_ice*gravity*(surface-z_g);
     1617                pressure       = ice_pressure + water_pressure;
     1618
     1619                for (int i=0;i<numnodes;i++){
     1620                        pe->values[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i];
     1621                        pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i];
     1622                }
     1623        }
     1624
     1625        /*Transform coordinate system*/
     1626        element->TransformLoadVectorCoord(pe,XYEnum);
     1627
     1628        /*Clean up and return*/
     1629        xDelete<IssmDouble>(basis);
     1630        xDelete<IssmDouble>(xyz_list);
     1631        xDelete<IssmDouble>(xyz_list_front);
     1632        delete gauss;
     1633        return pe;
    15731634}/*}}}*/
    15741635void StressbalanceAnalysis::GetBHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    … …  
    19371998                element->NodalFunctionsVelocity(vbasis,gauss);
    19381999
    1939                 element->NormalBase(&normal[0],xyz_list_base);
     2000                element->NormalSection(&normal[0],xyz_list_base);
    19402001                _assert_(normal[dim-1]<0.);
    19412002                bed_input->GetInputValue(&bed, gauss);

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

    r16838 r16839  
    6262                virtual void   NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss)=0;
    6363                virtual void   NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss)=0;
    64                 virtual void   NormalBase(IssmDouble* normal,IssmDouble* xyz_list)=0;
     64                virtual void   NormalSection(IssmDouble* normal,IssmDouble* xyz_list)=0;
    6565                virtual void   NormalTop(IssmDouble* normal,IssmDouble* xyz_list)=0;
    6666                virtual IssmDouble PureIceEnthalpy(IssmDouble pressure)=0;
    … …  
    111111                virtual void   GetVerticesCoordinatesBase(IssmDouble** xyz_list)=0;
    112112                virtual void   GetVerticesCoordinatesTop(IssmDouble** xyz_list)=0;
     113                virtual IssmDouble GetZcoord(Gauss* gauss)=0;
    113114                virtual void   GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0;
    114115
    … …  
    132133                virtual Gauss* NewGauss(int order)=0;
    133134      virtual Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order)=0;
     135      virtual Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert)=0;
    134136                virtual Gauss* NewGaussBase(int order)=0;
    135137                virtual Gauss* NewGaussTop(int order)=0;

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

    r16838 r16839  
    15111511}
    15121512/*}}}*/
    1513 /*FUNCTION Penta::GetQuadNormal {{{*/
    1514 void Penta:: GetQuadNormal(IssmDouble* normal,IssmDouble* xyz_list){
     1513/*FUNCTION Penta::NormalSection{{{*/
     1514void Penta::NormalSection(IssmDouble* normal,IssmDouble* xyz_list){
    15151515
    15161516        /*Build unit outward pointing vector*/
    … …  
    16341634/*}}}*/
    16351635/*FUNCTION Penta::GetZcoord {{{*/
    1636 IssmDouble Penta::GetZcoord(GaussPenta* gauss){
     1636IssmDouble Penta::GetZcoord(Gauss* gauss){
    16371637
    16381638        int    i;
    … …  
    24662466}
    24672467/*}}}*/
     2468/*FUNCTION Penta::JacobianDeterminantSurface{{{*/
     2469void Penta::JacobianDeterminantSurface(IssmDouble* pJdet,IssmDouble* xyz_list_quad,Gauss* gauss){
     2470
     2471        _assert_(gauss->Enum()==GaussPentaEnum);
     2472        this->GetQuadJacobianDeterminant(pJdet,xyz_list_quad,(GaussPenta*)gauss);
     2473
     2474}
     2475/*}}}*/
    24682476/*FUNCTION Penta::NoIceInElement {{{*/
    24692477bool   Penta::NoIceInElement(){
    … …  
    25422550Gauss* Penta::NewGauss(int order){
    25432551        return new GaussPenta(order,order);
     2552}
     2553/*}}}*/
     2554/*FUNCTION Penta::NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){{{*/
     2555Gauss* Penta::NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){
     2556
     2557        IssmDouble  area_coordinates[4][3];
     2558
     2559        GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,4);
     2560
     2561        return new GaussPenta(area_coordinates,order_horiz,order_vert);
    25442562}
    25452563/*}}}*/
    … …  
    31233141                return S;
    31243142        }
    3125 }
    3126 /*}}}*/
    3127 /*FUNCTION Penta::SurfaceNormal {{{*/
    3128 void Penta::SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]){
    3129 
    3130         int    i;
    3131         IssmDouble v13[3],v23[3];
    3132         IssmDouble normal[3];
    3133         IssmDouble normal_norm;
    3134 
    3135         for (i=0;i<3;i++){
    3136                 v13[i]=xyz_list[0][i]-xyz_list[2][i];
    3137                 v23[i]=xyz_list[1][i]-xyz_list[2][i];
    3138         }
    3139 
    3140         normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
    3141         normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
    3142         normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
    3143 
    3144         normal_norm=sqrt( pow(normal[0],2)+pow(normal[1],2)+pow(normal[2],2) );
    3145 
    3146         *(surface_normal)=normal[0]/normal_norm;
    3147         *(surface_normal+1)=normal[1]/normal_norm;
    3148         *(surface_normal+2)=normal[2]/normal_norm;
    31493143}
    31503144/*}}}*/
    … …  
    60726066
    60736067                /*Get normal vector to the bed */
    6074                 SurfaceNormal(&surface_normal[0],xyz_list_tria);
     6068                NormalTop(&surface_normal[0],&xyz_list_tria[0][0]);
    60756069
    60766070                bed_normal[0]=-surface_normal[0]; //Function is for upper surface, so the normal to the bed is the opposite of the result
    … …  
    81238117        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    81248118        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i+3][j];
    8125         SurfaceNormal(&surface_normal[0],xyz_list_tria);
     8119        NormalTop(&surface_normal[0],&xyz_list_tria[0][0]);
    81268120
    81278121        /* Start  looping on the number of gaussian points: */
    … …  
    87888782        ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum);
    87898783        GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],4);
    8790         GetQuadNormal(&normal[0],xyz_list_front);
     8784        NormalSection(&normal[0],xyz_list_front);
    87918785
    87928786        /*Initialize gauss points*/
    … …  
    88218815        /*Clean up and return*/
    88228816        xDelete<IssmDouble>(basis);
     8817        xDelete<IssmDouble>(xyz_list_front);
    88238818        delete gauss;
    88248819        return pe;
    … …  
    88838878        ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum);
    88848879        GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],4);
    8885         GetQuadNormal(&normal[0],xyz_list_front);
     8880        NormalSection(&normal[0],xyz_list_front);
    88868881
    88878882        /*Start looping on Gaussian points*/

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

    r16838 r16839  
    9797                IssmDouble GetMaterialParameter(int enum_in);
    9898                void   GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type);
    99                 IssmDouble GetZcoord(GaussPenta* gauss);
     99                IssmDouble GetZcoord(Gauss* gauss);
    100100                void   GetVectorFromInputs(Vector<IssmDouble>* vector,int name_enum);
    101101                void   GetVerticesCoordinates(IssmDouble** pxyz_list);
    … …  
    199199                void           AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum);
    200200                void             NormalBase(IssmDouble* bed_normal, IssmDouble* xyz_list);
     201                void           NormalSection(IssmDouble* normal,IssmDouble* xyz_list);
    201202                void             NormalTop(IssmDouble* bed_normal, IssmDouble* xyz_list);
    202203                ElementMatrix* CreateBasalMassMatrix(void);
    … …  
    232233                void           GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype);
    233234                void             GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity);
    234                 void           GetQuadNormal(IssmDouble* normal,IssmDouble* xyz_list);
    235235                void           GetStrainRate3dHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss, Input* vx_input, Input* vy_input);
    236236                void           GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, Gauss* gauss, Input* vx_input, Input* vy_input, Input* vz_input);
    … …  
    247247                bool           IsNodeOnShelfFromFlags(IssmDouble* flags);
    248248                void           JacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
    249                 void           JacobianDeterminantSurface(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
     249                void           JacobianDeterminantSurface(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss);
    250250                void           JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
    251251                void           JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
    … …  
    254254                Gauss*         NewGauss(int order);
    255255      Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){_error_("not implemented yet");};
     256      Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert);
    256257                Gauss*         NewGaussBase(int order);
    257258                Gauss*         NewGaussTop(int order);
    … …  
    266267                void             SetClone(int* minranks);
    267268                Tria*            SpawnTria(int location);
    268                 void             SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]);
    269269                IssmDouble     StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa);
    270270                void           TransformLoadVectorCoord(ElementVector* pe,int transformenum);

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

    r16838 r16839  
    124124                void        NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    125125                bool        NoIceInElement(){_error_("not implemented yet");};
    126                 void        NormalBase(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
     126                void        NormalSection(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
    127127                void        NormalTop(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
    128128                int         NumberofNodesVelocity(void){_error_("not implemented yet");};
    … …  
    141141                void        GetInputValue(IssmDouble* pvalue,int enum_type){_error_("not implemented yet");};
    142142                void        GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype){_error_("not implemented yet");};
     143                IssmDouble  GetZcoord(Gauss* gauss){_error_("not implemented yet");};
    143144                Gauss*      NewGauss(void){_error_("not implemented yet");};
    144145                Gauss*      NewGauss(int order){_error_("not implemented yet");};
    145146      Gauss*      NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){_error_("not implemented yet");};
     147      Gauss*      NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){_error_("not implemented yet");};
    146148                Gauss*      NewGaussBase(int order){_error_("not implemented yet");};
    147149                Gauss*      NewGaussTop(int order){_error_("not implemented yet");};

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

    r16838 r16839  
    11681168
    11691169}/*}}}*/
    1170 /*FUNCTION Tria::NormalBase {{{*/
    1171 void Tria::NormalBase(IssmDouble* normal,IssmDouble* xyz_list){
     1170/*FUNCTION Tria::NormalSection{{{*/
     1171void Tria::NormalSection(IssmDouble* normal,IssmDouble* xyz_list){
    11721172
    11731173        /*Build unit outward pointing vector*/
    … …  
    26972697}
    26982698/*}}}*/
    2699 /*FUNCTION Tria::SurfaceNormal{{{*/
    2700 void Tria::SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]){
    2701 
    2702         IssmDouble v13[3],v23[3];
    2703         IssmDouble normal[3];
    2704         IssmDouble normal_norm;
    2705 
    2706         for(int i=0;i<3;i++){
    2707                 v13[i]=xyz_list[0][i]-xyz_list[2][i];
    2708                 v23[i]=xyz_list[1][i]-xyz_list[2][i];
    2709         }
    2710 
    2711         normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
    2712         normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
    2713         normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
    2714 
    2715         normal_norm=sqrt( normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
    2716 
    2717         *(surface_normal+0) = normal[0]/normal_norm;
    2718         *(surface_normal+1) = normal[1]/normal_norm;
    2719         *(surface_normal+2) = normal[2]/normal_norm;
    2720 }
    2721 /*}}}*/
    27222699/*FUNCTION Tria::TimeAdapt{{{*/
    27232700IssmDouble  Tria::TimeAdapt(void){
    … …  
    39763953        ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum);
    39773954        GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],2);
    3978         NormalBase(&normal[0],xyz_list_front);
     3955        NormalSection(&normal[0],xyz_list_front);
    39793956
    39803957        /*Start looping on Gaussian points*/
    … …  
    41224099                GetNodalFunctionsVelocity(vbasis, gauss);
    41234100
    4124                 NormalBase(&normal[0],&xyz_list_seg[0][0]);
     4101                NormalSection(&normal[0],&xyz_list_seg[0][0]);
    41254102                _assert_(normal[1]<0.);
    41264103                bed_input->GetInputValue(&bed, gauss);
    … …  
    42424219        ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum);
    42434220        GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],2);
    4244         NormalBase(&normal[0],xyz_list_front);
     4221        NormalSection(&normal[0],xyz_list_front);
    42454222
    42464223        /*Start looping on Gaussian points*/
    … …  
    72207197        this->EdgeOnSurfaceIndices(&indices[0],&indices[1]);
    72217198        for(int i=0;i<NUMVERTICES1D;i++) for(int j=0;j<2;j++) xyz_list_seg[i][j]=xyz_list[indices[i]][j];
    7222         NormalBase(&normal[0],&xyz_list_seg[0][0]);
     7199        NormalSection(&normal[0],&xyz_list_seg[0][0]);
    72237200
    72247201        /* Start  looping on the number of gaussian points: */
    … …  
    72707247        this->EdgeOnBedIndices(&indices[0],&indices[1]);
    72717248        for(int i=0;i<NUMVERTICES1D;i++) for(int j=0;j<2;j++) xyz_list_seg[i][j]=xyz_list[indices[i]][j];
    7272         NormalBase(&normal[0],&xyz_list_seg[0][0]);
     7249        NormalSection(&normal[0],&xyz_list_seg[0][0]);
    72737250
    72747251        /* Start  looping on the number of gaussian points: */

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

    r16838 r16839  
    257257                void           GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
    258258                IssmDouble     GetGroundedPortion(IssmDouble* xyz_list);
    259                 void           NormalBase(IssmDouble* normal,IssmDouble* xyz_list);
     259                IssmDouble     GetZcoord(Gauss* gauss){_error_("not implemented");};
     260                void           NormalSection(IssmDouble* normal,IssmDouble* xyz_list);
    260261                void           NormalTop(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
    261262                IssmDouble     GetMaterialParameter(int enum_in);
    … …  
    284285                Gauss*         NewGauss(int order);
    285286      Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order);
     287      Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){_error_("not implemented yet");};
    286288                Gauss*         NewGaussBase(int order){_error_("not implemented yet");};
    287289                Gauss*         NewGaussTop(int order){_error_("not implemented yet");};
    … …  
    296298                Seg*             SpawnSeg(int index1,int index2);
    297299                IssmDouble     StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");};
    298                 void             SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]);
    299300                void           TransformLoadVectorCoord(ElementVector* pe,int transformenum);
    300301                void           TransformLoadVectorCoord(ElementVector* pe,int* transformenum_list);
Note: See TracChangeset for help on using the changeset viewer.