Changeset 16844


Ignore:
Timestamp:
11/21/13 07:30:30 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: FS front

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

Legend:

Unmodified
Added
Removed

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

    r16839 r16844  
    11141114
    11151115        /*Intermediaries*/
    1116         IssmDouble  rho_ice,rho_water,gravity;
    11171116        IssmDouble  Jdet,thickness,bed,water_pressure,ice_pressure;
    11181117        IssmDouble  surface_under_water,base_under_water,pressure;
    … …  
    11291128
    11301129        /*Retrieve all inputs and parameters*/
    1131         Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input);
    1132         Input* bed_input      =element->GetInput(BedEnum);       _assert_(bed_input);
    1133         rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum);
    1134         rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
    1135         gravity   = element->GetMaterialParameter(ConstantsGEnum);
     1130        Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
     1131        Input* bed_input       = element->GetInput(BedEnum);       _assert_(bed_input);
     1132        IssmDouble rho_water  = element->GetMaterialParameter(MaterialsRhoWaterEnum);
     1133        IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
     1134        IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
    11361135        element->GetVerticesCoordinates(&xyz_list);
    11371136        element->ZeroLevelsetCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
    … …  
    15741573
    15751574        /*Intermediaries*/
    1576         IssmDouble  rho_ice,rho_water,gravity;
    1577         IssmDouble  Jdet,surface,z_g,water_pressure,ice_pressure;
     1575        IssmDouble  Jdet,surface,z,water_pressure,ice_pressure;
    15781576        IssmDouble  surface_under_water,base_under_water,pressure;
    15791577        IssmDouble* xyz_list       = NULL;
    … …  
    15831581
    15841582        /*Fetch number of nodes and dof for this finite element*/
    1585         int numnodes = element->GetNumberOfNodes();
     1583        int numnodes    = element->GetNumberOfNodes();
     1584        int numvertices = element->GetNumberOfVertices();
    15861585
    15871586        /*Initialize Element vector and other vectors*/
    … …  
    15901589
    15911590        /*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);
     1591        Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input);
     1592        IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum);
     1593        IssmDouble rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
     1594        IssmDouble gravity   = element->GetMaterialParameter(ConstantsGEnum);
    15961595        element->GetVerticesCoordinates(&xyz_list);
    15971596        element->ZeroLevelsetCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
    … …  
    15991598
    16001599        /*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];
     1600        IssmDouble zmax=xyz_list[0*3+2]; for(int i=1;i<numvertices;i++) if(xyz_list[i*3+2]>zmax) zmax=xyz_list[i*3+2];
     1601        IssmDouble zmin=xyz_list[0*3+2]; for(int i=1;i<numvertices;i++) if(xyz_list[i*3+2]<zmin) zmin=xyz_list[i*3+2];
    16031602        if(zmax>0 && zmin<0) gauss=element->NewGauss(xyz_list,xyz_list_front,3,10);//refined in vertical because of the sea level discontinuity
    16041603        else                 gauss=element->NewGauss(xyz_list,xyz_list_front,3,3);
    … …  
    16091608                gauss->GaussPoint(ig);
    16101609                surface_input->GetInputValue(&surface,gauss);
    1611                 z_g=element->GetZcoord(gauss);
     1610                z=element->GetZcoord(gauss);
    16121611                element->NodalFunctions(basis,gauss);
    16131612                element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
    16141613
    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);
     1614                water_pressure = rho_water*gravity*min(0.,z);//0 if the gaussian point is above water level
     1615                ice_pressure   = rho_ice*gravity*(surface-z);
    16171616                pressure       = ice_pressure + water_pressure;
    16181617
    … …  
    19981997                element->NodalFunctionsVelocity(vbasis,gauss);
    19991998
    2000                 element->NormalSection(&normal[0],xyz_list_base);
     1999                element->NormalBase(&normal[0],xyz_list_base);
    20012000                _assert_(normal[dim-1]<0.);
    20022001                bed_input->GetInputValue(&bed, gauss);
    … …  
    20242023ElementVector* StressbalanceAnalysis::CreatePVectorFSFront(Element* element){/*{{{*/
    20252024
    2026         return NULL;
     2025        /*If no front, return NULL*/
     2026        if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) return NULL;
     2027
     2028        /*Intermediaries*/
     2029        int         i,meshtype,dim;
     2030        IssmDouble  Jdet,pressure,surface,z;
     2031        IssmDouble      normal[3];
     2032        IssmDouble *xyz_list       = NULL;
     2033        IssmDouble *xyz_list_front = NULL;
     2034        Gauss      *gauss          = NULL;
     2035
     2036        /*Get basal element*/
     2037        if(!element->IsOnBed() || !element->IsFloating()) return NULL;
     2038
     2039        /*Get problem dimension*/
     2040        element->FindParam(&meshtype,MeshTypeEnum);
     2041        switch(meshtype){
     2042                case Mesh2DverticalEnum: dim = 2; break;
     2043                case Mesh3DEnum:         dim = 3; break;
     2044                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     2045        }
     2046
     2047        /*Fetch number of nodes and dof for this finite element*/
     2048        int vnumnodes   = element->GetNumberOfNodesVelocity();
     2049        int pnumnodes   = element->GetNumberOfNodesPressure();
     2050        int numvertices = element->GetNumberOfVertices();
     2051
     2052        /*Prepare coordinate system list*/
     2053        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     2054        if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
     2055        else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
     2056        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     2057
     2058        /*Initialize vectors*/
     2059        ElementVector* pe     = element->NewElementVector(FSvelocityEnum);
     2060        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
     2061
     2062        /*Retrieve all inputs and parameters*/
     2063        element->GetVerticesCoordinates(&xyz_list);
     2064        element->ZeroLevelsetCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
     2065        element->NormalSection(&normal[0],xyz_list_front);
     2066        Input* surface_input  = element->GetInput(SurfaceEnum); _assert_(surface_input);
     2067        IssmDouble  rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum);
     2068        IssmDouble  gravity   = element->GetMaterialParameter(ConstantsGEnum);
     2069
     2070        /*Initialize gauss points*/
     2071        IssmDouble zmax=xyz_list[0*3+2]; for(int i=1;i<numvertices;i++) if(xyz_list[i*3+2]>zmax) zmax=xyz_list[i*3+2];
     2072        IssmDouble zmin=xyz_list[0*3+2]; for(int i=1;i<numvertices;i++) if(xyz_list[i*3+2]<zmin) zmin=xyz_list[i*3+2];
     2073        if(zmax>0 && zmin<0) gauss=element->NewGauss(xyz_list,xyz_list_front,3,30);//refined in vertical because of the sea level discontinuity
     2074        else                 gauss=element->NewGauss(xyz_list,xyz_list_front,3,3);
     2075
     2076        /* Start  looping on the number of gaussian points: */
     2077        for(int ig=gauss->begin();ig<gauss->end();ig++){
     2078                gauss->GaussPoint(ig);
     2079
     2080                element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
     2081                element->NodalFunctionsVelocity(vbasis,gauss);
     2082                surface_input->GetInputValue(&surface,gauss);
     2083                z=element->GetZcoord(gauss);
     2084                pressure = rho_water*gravity*min(0.,z);//0 if the gaussian point is above water level
     2085
     2086                for (int i=0;i<vnumnodes;i++){
     2087                        pe->values[dim*i+0]+= pressure*Jdet*gauss->weight*normal[0]*vbasis[i];
     2088                        pe->values[dim*i+1]+= pressure*Jdet*gauss->weight*normal[1]*vbasis[i];
     2089                        if(dim==3) pe->values[dim*i+2]+= pressure*Jdet*gauss->weight*normal[2]*vbasis[i];
     2090                }
     2091        }
     2092
     2093        /*Transform coordinate system*/
     2094        element->TransformLoadVectorCoord(pe,cs_list);
     2095
     2096        /*Clean up and return*/
     2097        delete gauss;
     2098        xDelete<int>(cs_list);
     2099        xDelete<IssmDouble>(vbasis);
     2100        xDelete<IssmDouble>(xyz_list);
     2101        xDelete<IssmDouble>(xyz_list_front);
     2102        return pe;
    20272103}/*}}}*/
    20282104void StressbalanceAnalysis::GetBFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/

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

    r16840 r16844  
    6565                virtual void   NormalSection(IssmDouble* normal,IssmDouble* xyz_list)=0;
    6666                virtual void   NormalTop(IssmDouble* normal,IssmDouble* xyz_list)=0;
     67                virtual void   NormalBase(IssmDouble* normal,IssmDouble* xyz_list)=0;
    6768                virtual IssmDouble PureIceEnthalpy(IssmDouble pressure)=0;
    6869                virtual void   TransformLoadVectorCoord(ElementVector* pe,int transformenum)=0;

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

    r16840 r16844  
    127127                void        NormalSection(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
    128128                void        NormalTop(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
     129                void        NormalBase(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
    129130                int         NumberofNodesVelocity(void){_error_("not implemented yet");};
    130131                int         NumberofNodesPressure(void){_error_("not implemented yet");};

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

    r16840 r16844  
    261261                void           NormalSection(IssmDouble* normal,IssmDouble* xyz_list);
    262262                void           NormalTop(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
     263                void           NormalBase(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
    263264                IssmDouble     GetMaterialParameter(int enum_in);
    264265                Input*         GetInput(int inputenum);
Note: See TracChangeset for help on using the changeset viewer.