Changeset 16807


Ignore:
Timestamp:
11/16/13 20:14:26 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: done with FS shelf

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

Legend:

Unmodified
Added
Removed

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

    r16805 r16807  
    431431                        iomodel->FetchData(&z,NULL,NULL,MeshZEnum);
    432432                        switch(finiteelement){
     433                                case P1Enum:
     434                                        for(i=0;i<iomodel->numberofvertices;i++){
     435                                                if(iomodel->my_vertices[i]){
     436                                                        if(reCast<int,IssmDouble>(vertices_type[i])==NoneApproximationEnum){
     437                                                                constraints->AddObject(new SpcStatic(count+1,iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+i+1,1,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum));
     438                                                                count++;
     439                                                        }
     440                                                }
     441                                        }
     442                                        break;
    433443                                case P1bubbleEnum:
    434444                                        for(i=0;i<iomodel->numberofvertices;i++){
    … …  
    909919ElementVector* StressbalanceAnalysis::CreatePVectorFSShelf(Element* element){/*{{{*/
    910920
    911         return NULL;
     921        int         i,meshtype,dim;
     922        IssmDouble  Jdet,water_pressure,bed;
     923        IssmDouble      normal[3];
     924        IssmDouble *xyz_list_base = NULL;
     925
     926        /*Get basal element*/
     927        if(!element->IsOnBed() || !element->IsFloating()) return NULL;
     928
     929        /*Get problem dimension*/
     930        element->FindParam(&meshtype,MeshTypeEnum);
     931        switch(meshtype){
     932                case Mesh2DverticalEnum: dim = 2; break;
     933                case Mesh3DEnum:         dim = 3; break;
     934                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     935        }
     936
     937        /*Fetch number of nodes and dof for this finite element*/
     938        int vnumnodes = element->GetNumberOfNodesVelocity();
     939        int pnumnodes = element->GetNumberOfNodesPressure();
     940
     941        /*Prepare coordinate system list*/
     942        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     943        if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
     944        else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
     945        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     946
     947        /*Initialize vectors*/
     948        ElementVector* pe     = element->NewElementVector(FSvelocityEnum);
     949        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
     950
     951        /*Retrieve all inputs and parameters*/
     952        element->GetVerticesCoordinatesBase(&xyz_list_base);
     953        Input*      bed_input=element->GetInput(BedEnum); _assert_(bed_input);
     954        IssmDouble  rho_water=element->GetMaterialParameter(MaterialsRhoWaterEnum);
     955        IssmDouble  gravity  =element->GetMaterialParameter(ConstantsGEnum);
     956
     957        /* Start  looping on the number of gaussian points: */
     958        Gauss* gauss=element->NewGaussBase(5);
     959        for(int ig=gauss->begin();ig<gauss->end();ig++){
     960                gauss->GaussPoint(ig);
     961
     962                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
     963                element->NodalFunctionsVelocity(vbasis,gauss);
     964
     965                element->NormalBase(&normal[0],xyz_list_base);
     966                _assert_(normal[dim-1]<0.);
     967                bed_input->GetInputValue(&bed, gauss);
     968                water_pressure=gravity*rho_water*bed;
     969
     970                for(i=0;i<vnumnodes;i++){
     971                        pe->values[i*dim+0] += water_pressure*gauss->weight*Jdet*vbasis[i]*normal[0];
     972                        pe->values[i*dim+1] += water_pressure*gauss->weight*Jdet*vbasis[i]*normal[1];
     973                        if(dim==3){
     974                                pe->values[i*dim+2]+=water_pressure*gauss->weight*Jdet*vbasis[i]*normal[2];
     975                        }
     976                }
     977        }
     978
     979        /*Transform coordinate system*/
     980        element->TransformLoadVectorCoord(pe,cs_list);
     981
     982        /*Clean up and return*/
     983        delete gauss;
     984        xDelete<int>(cs_list);
     985        xDelete<IssmDouble>(vbasis);
     986        xDelete<IssmDouble>(xyz_list_base);
     987        return pe;
    912988}/*}}}*/
    913989ElementVector* StressbalanceAnalysis::CreatePVectorFSFront(Element* element){/*{{{*/

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

    r16805 r16807  
    5656                virtual void   NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss)=0;
    5757                virtual void   NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss)=0;
     58                virtual void   NormalBase(IssmDouble* normal,IssmDouble* xyz_list)=0;
    5859                virtual void   TransformLoadVectorCoord(ElementVector* pe,int transformenum)=0;
    5960                virtual void   TransformLoadVectorCoord(ElementVector* pe,int* transformenum_list)=0;

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

    r16805 r16807  
    162162/*}}}*/
    163163
    164 /*FUNCTION Penta::BedNormal {{{*/
    165 void Penta::BedNormal(IssmDouble* bed_normal, IssmDouble xyz_list[3][3]){
     164/*FUNCTION Penta::NormalBase {{{*/
     165void Penta::NormalBase(IssmDouble* bed_normal,IssmDouble* xyz_list){
    166166
    167167        int i;
    … …  
    171171
    172172        for (i=0;i<3;i++){
    173                 v13[i]=xyz_list[0][i]-xyz_list[2][i];
    174                 v23[i]=xyz_list[1][i]-xyz_list[2][i];
     173                v13[i]=xyz_list[0*3+i]-xyz_list[2*3+i];
     174                v23[i]=xyz_list[1*3+i]-xyz_list[2*3+i];
    175175        }
    176176
    … …  
    181181
    182182        /*Bed normal is opposite to surface normal*/
    183         *(bed_normal)=-normal[0]/normal_norm;
    184         *(bed_normal+1)=-normal[1]/normal_norm;
    185         *(bed_normal+2)=-normal[2]/normal_norm;
     183        bed_normal[0]=-normal[0]/normal_norm;
     184        bed_normal[1]=-normal[1]/normal_norm;
     185        bed_normal[2]=-normal[2]/normal_norm;
    186186}
    187187/*}}}*/
    … …  
    304304
    305305                /*Get normal vector to the bed */
    306                 BedNormal(&bed_normal[0],xyz_list_tria);
     306                NormalBase(&bed_normal[0],&xyz_list_tria[0][0]);
    307307
    308308                /*basalforce*/
    … …  
    47884788
    47894789                        /*geothermal heatflux*/
    4790                         BedNormal(&normal_base[0],xyz_list_tria);
     4790                        NormalBase(&normal_base[0],&xyz_list_tria[0][0]);
    47914791                        heatflux=0.;
    47924792                        for(i=0;i<3;i++) heatflux+=(vec_heatflux[i])*normal_base[i];
    … …  
    68886888                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    68896889
    6890                 BedNormal(&bed_normal[0],xyz_list_tria);
     6890                NormalBase(&bed_normal[0],&xyz_list_tria[0][0]);
    68916891                friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum);
    68926892
    … …  
    80708070                vzSSA_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
    80718071
    8072                 BedNormal(&bed_normal[0],xyz_list_tria);
     8072                NormalBase(&bed_normal[0],&xyz_list_tria[0][0]);
    80738073                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    80748074                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    … …  
    82358235                vzHO_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
    82368236
    8237                 BedNormal(&bed_normal[0],xyz_list_tria);
     8237                NormalBase(&bed_normal[0],&xyz_list_tria[0][0]);
    82388238                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    82398239                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    … …  
    89348934                GetNodalFunctionsVelocity(vbasis, gauss);
    89358935
    8936                 BedNormal(&bed_normal[0],xyz_list_tria);
     8936                NormalBase(&bed_normal[0],&xyz_list_tria[0][0]);
    89378937                bed_input->GetInputValue(&bed, gauss);
    89388938                if(shelf_dampening){ //add dampening to avoid too high vertical velocities when not in hydrostatic equilibrium

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

    r16805 r16807  
    192192                void           AddInput(int input_enum, IssmDouble* values, int interpolation_enum);
    193193                void           AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum);
    194                 void             BedNormal(IssmDouble* bed_normal, IssmDouble xyz_list[3][3]);
     194                void             NormalBase(IssmDouble* bed_normal, IssmDouble* xyz_list);
    195195                ElementMatrix* CreateBasalMassMatrix(void);
    196196                ElementMatrix* CreateKMatrix(void);

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

    r16805 r16807  
    114114                void        NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");};
    115115                bool        NoIceInElement(){_error_("not implemented yet");};
     116                void        NormalBase(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
    116117                int         NumberofNodesVelocity(void){_error_("not implemented yet");};
    117118                int         NumberofNodesPressure(void){_error_("not implemented yet");};

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

    r16805 r16807  
    11621162
    11631163}/*}}}*/
    1164 /*FUNCTION Tria::GetSegmentNormal {{{*/
    1165 void Tria:: GetSegmentNormal(IssmDouble* normal,IssmDouble xyz_list[2][3]){
     1164/*FUNCTION Tria::NormalBase {{{*/
     1165void Tria::NormalBase(IssmDouble* normal,IssmDouble* xyz_list){
    11661166
    11671167        /*Build unit outward pointing vector*/
    … …  
    11691169        IssmDouble norm;
    11701170
    1171         vector[0]=xyz_list[1][0] - xyz_list[0][0];
    1172         vector[1]=xyz_list[1][1] - xyz_list[0][1];
     1171        vector[0]=xyz_list[1*3+0] - xyz_list[0*3+0];
     1172        vector[1]=xyz_list[1*3+1] - xyz_list[0*3+1];
    11731173
    11741174        norm=sqrt(vector[0]*vector[0] + vector[1]*vector[1]);
    … …  
    19761976
    19771977        bool onbed;
    1978         inputs->GetInputValue(&onbed,MeshElementonbedEnum);
     1978        int meshtype;
     1979        this->parameters->FindParam(&meshtype,MeshTypeEnum);
     1980        switch(meshtype){
     1981                case Mesh2DverticalEnum:
     1982                        return HasEdgeOnBed();
     1983                        break;
     1984                case Mesh2DhorizontalEnum:
     1985                        inputs->GetInputValue(&onbed,MeshElementonbedEnum);
     1986                        break;
     1987                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     1988        }
    19791989        return onbed;
    19801990}
    … …  
    38763886        GetZeroLevelsetCoordinates(&xyz_list_front[0][0],xyz_list,MaskIceLevelsetEnum);
    38773887        GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,2);
    3878         GetSegmentNormal(&normal[0],xyz_list_front);
     3888        NormalBase(&normal[0],&xyz_list_front[0][0]);
    38793889
    38803890        /*Start looping on Gaussian points*/
    … …  
    40224032                GetNodalFunctionsVelocity(vbasis, gauss);
    40234033
    4024                 GetSegmentNormal(&normal[0],xyz_list_seg);
     4034                NormalBase(&normal[0],&xyz_list_seg[0][0]);
    40254035                _assert_(normal[1]<0.);
    40264036                bed_input->GetInputValue(&bed, gauss);
    … …  
    41574167        GetZeroLevelsetCoordinates(&xyz_list_front[0][0],xyz_list,MaskIceLevelsetEnum);
    41584168        GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,2);
    4159         GetSegmentNormal(&normal[0],xyz_list_front);
     4169        NormalBase(&normal[0],&xyz_list_front[0][0]);
    41604170
    41614171        /*Start looping on Gaussian points*/
    … …  
    71257135        this->EdgeOnSurfaceIndices(&indices[0],&indices[1]);
    71267136        for(int i=0;i<NUMVERTICES1D;i++) for(int j=0;j<2;j++) xyz_list_seg[i][j]=xyz_list[indices[i]][j];
    7127         GetSegmentNormal(&normal[0],xyz_list_seg);
     7137        NormalBase(&normal[0],&xyz_list_seg[0][0]);
    71287138
    71297139        /* Start  looping on the number of gaussian points: */
    … …  
    71757185        this->EdgeOnBedIndices(&indices[0],&indices[1]);
    71767186        for(int i=0;i<NUMVERTICES1D;i++) for(int j=0;j<2;j++) xyz_list_seg[i][j]=xyz_list[indices[i]][j];
    7177         GetSegmentNormal(&normal[0],xyz_list_seg);
     7187        NormalBase(&normal[0],&xyz_list_seg[0][0]);
    71787188
    71797189        /* Start  looping on the number of gaussian points: */

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

    r16805 r16807  
    248248                void           GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
    249249                IssmDouble     GetGroundedPortion(IssmDouble* xyz_list);
    250                 void           GetSegmentNormal(IssmDouble* normal,IssmDouble xyz_list[2][3]);
     250                void           NormalBase(IssmDouble* normal,IssmDouble* xyz_list);
    251251                IssmDouble     GetMaterialParameter(int enum_in);
    252252                void           GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[3][3],int levelsetenum);
Note: See TracChangeset for help on using the changeset viewer.