Changeset 16010


Ignore:
Timestamp:
08/29/13 09:41:24 (12 years ago)
Author:
Mathieu Morlighem
Message:

NEW: added smoothed surface, required for bal vel

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

Legend:

Unmodified
Added
Removed

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

    r16007 r16010  
    2020
    2121        if(VerboseSolution()) _printf0_("call computational core:\n");
    22         surfaceslope_core(femmodel);
     22        femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SmoothedSurfaceSlopeXAnalysisEnum);
     23        solutionsequence_linear(femmodel);
     24        femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SmoothedSurfaceSlopeYAnalysisEnum);
     25        solutionsequence_linear(femmodel);
    2326
    2427        if(VerboseSolution()) _printf0_("call computational core:\n");
    … …  
    2831        if(save_results){
    2932                if(VerboseSolution()) _printf0_("   saving results\n");
     33                InputToResultx(femmodel,SurfaceSlopeXEnum);
     34                InputToResultx(femmodel,SurfaceSlopeYEnum);
    3035                InputToResultx(femmodel,VelEnum);
    3136        }

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

    r16007 r16010  
    247247                        return CreateKMatrixBalancevelocity();
    248248                        break;
     249                case SmoothedSurfaceSlopeXAnalysisEnum: case SmoothedSurfaceSlopeYAnalysisEnum:
     250                        return CreateKMatrixSmoothedSlope();
     251                        break;
    249252                #endif
    250253                #ifdef _HAVE_CONTROL_
    … …  
    391394                case BalancevelocityAnalysisEnum:
    392395                        return CreatePVectorBalancevelocity();
     396                        break;
     397                case SmoothedSurfaceSlopeXAnalysisEnum:
     398                        return CreatePVectorSmoothedSlopeX();
     399                        break;
     400                case SmoothedSurfaceSlopeYAnalysisEnum:
     401                        return CreatePVectorSmoothedSlopeY();
    393402                        break;
    394403#endif
    … …  
    16081617                case BalancevelocityAnalysisEnum:
    16091618                        InputUpdateFromSolutionOneDof(solution,VelEnum);
     1619                        break;
     1620                case SmoothedSurfaceSlopeXAnalysisEnum:
     1621                        InputUpdateFromSolutionOneDof(solution,SurfaceSlopeXEnum);
     1622                        break;
     1623                case SmoothedSurfaceSlopeYAnalysisEnum:
     1624                        InputUpdateFromSolutionOneDof(solution,SurfaceSlopeYEnum);
    16101625                        break;
    16111626                #endif
    … …  
    72967311}
    72977312/*}}}*/
     7313/*FUNCTION Tria::CreateKMatrixSmoothedSlope {{{*/
     7314ElementMatrix* Tria::CreateKMatrixSmoothedSlope(void){
     7315
     7316        /* Intermediaries */
     7317        IssmDouble D_scalar,Jdet,thickness;
     7318        IssmDouble xyz_list[NUMVERTICES][3];
     7319        IssmDouble D[2][2];
     7320        IssmDouble l=12.;
     7321
     7322        /*Fetch number of nodes and dof for this finite element*/
     7323        int numnodes = this->NumberofNodes();
     7324
     7325        /*Initialize Element matrix and vectors*/
     7326        ElementMatrix* Ke    = CreateMassMatrix();
     7327        IssmDouble*    B     = xNew<IssmDouble>(2*numnodes);
     7328
     7329        /*Retrieve all inputs and parameters*/
     7330        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     7331        Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
     7332
     7333        /* Start looping on the number of gaussian points: */
     7334        GaussTria* gauss=new GaussTria(2);
     7335        for(int ig=gauss->begin();ig<gauss->end();ig++){
     7336
     7337                gauss->GaussPoint(ig);
     7338
     7339                GetBMasstransport(B,&xyz_list[0][0], gauss);
     7340                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     7341                thickness_input->GetInputValue(&thickness,gauss);
     7342
     7343                D_scalar=gauss->weight*Jdet*(l*thickness)*(l*thickness);
     7344                D[0][0]=D_scalar;
     7345           D[1][0]=0.;
     7346                D[0][1]=0.;
     7347                D[1][1]=D_scalar;
     7348
     7349                TripleMultiply(B,2,numnodes,1,
     7350                                        &D[0][0],2,2,0,
     7351                                        B,2,numnodes,0,
     7352                                        &Ke->values[0],1);
     7353        }
     7354
     7355        /*Clean up and return*/
     7356        delete gauss;
     7357        xDelete<IssmDouble>(B);
     7358        return Ke;
     7359}
     7360/*}}}*/
    72987361/*FUNCTION Tria::CreatePVectorBalancethickness{{{*/
    72997362ElementVector* Tria::CreatePVectorBalancethickness(void){
    … …  
    74777540        xDelete<IssmDouble>(HNx);
    74787541        xDelete<IssmDouble>(HNy);
     7542        delete gauss;
     7543        return pe;
     7544}
     7545/*}}}*/
     7546/*FUNCTION Tria::CreatePVectorSmoothedSlopeX{{{*/
     7547ElementVector* Tria::CreatePVectorSmoothedSlopeX(void){
     7548
     7549        /*Intermediaries */
     7550        IssmDouble xyz_list[NUMVERTICES][3];
     7551        IssmDouble Jdettria;
     7552        IssmDouble thickness,slope[2];
     7553        IssmDouble taud_x;
     7554
     7555        /*Fetch number of nodes and dof for this finite element*/
     7556        int numnodes = this->NumberofNodes();
     7557
     7558        /*Initialize Element vector and other vectors*/
     7559        ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
     7560        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     7561
     7562        /*Retrieve all inputs and parameters*/
     7563        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     7564        Input* H_input       = inputs->GetInput(ThicknessEnum); _assert_(H_input);
     7565        Input* surface_input = inputs->GetInput(SurfaceEnum);   _assert_(surface_input);
     7566
     7567        /* Start  looping on the number of gaussian points: */
     7568        GaussTria* gauss=new GaussTria(2);
     7569        for(int ig=gauss->begin();ig<gauss->end();ig++){
     7570
     7571                gauss->GaussPoint(ig);
     7572
     7573                H_input->GetInputValue(&thickness,gauss);
     7574                surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
     7575                taud_x = matpar->GetRhoIce()*matpar->GetG()*thickness*slope[0];
     7576
     7577                GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
     7578                GetNodalFunctions(basis,gauss);
     7579
     7580                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*taud_x*basis[i];
     7581        }
     7582
     7583        /*Clean up and return*/
     7584        xDelete<IssmDouble>(basis);
     7585        delete gauss;
     7586        return pe;
     7587}
     7588/*}}}*/
     7589/*FUNCTION Tria::CreatePVectorSmoothedSlopeY{{{*/
     7590ElementVector* Tria::CreatePVectorSmoothedSlopeY(void){
     7591
     7592        /*Intermediaries */
     7593        IssmDouble xyz_list[NUMVERTICES][3];
     7594        IssmDouble Jdettria;
     7595        IssmDouble thickness,slope[2];
     7596        IssmDouble taud_y;
     7597
     7598        /*Fetch number of nodes and dof for this finite element*/
     7599        int numnodes = this->NumberofNodes();
     7600
     7601        /*Initialize Element vector and other vectors*/
     7602        ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
     7603        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     7604
     7605        /*Retrieve all inputs and parameters*/
     7606        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     7607        Input* H_input       = inputs->GetInput(ThicknessEnum); _assert_(H_input);
     7608        Input* surface_input = inputs->GetInput(SurfaceEnum);   _assert_(surface_input);
     7609
     7610        /* Start  looping on the number of gaussian points: */
     7611        GaussTria* gauss=new GaussTria(2);
     7612        for(int ig=gauss->begin();ig<gauss->end();ig++){
     7613
     7614                gauss->GaussPoint(ig);
     7615
     7616                H_input->GetInputValue(&thickness,gauss);
     7617                surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
     7618                taud_y = matpar->GetRhoIce()*matpar->GetG()*thickness*slope[1];
     7619
     7620                GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
     7621                GetNodalFunctions(basis,gauss);
     7622
     7623                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*taud_y*basis[i];
     7624        }
     7625
     7626        /*Clean up and return*/
     7627        xDelete<IssmDouble>(basis);
    74797628        delete gauss;
    74807629        return pe;

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

    r16007 r16010  
    184184                ElementMatrix* CreateKMatrixBalancethickness_CG(void);
    185185                ElementMatrix* CreateKMatrixBalancevelocity(void);
     186                ElementMatrix* CreateKMatrixSmoothedSlope(void);
    186187                ElementMatrix* CreateKMatrixMelting(void);
    187188                ElementMatrix* CreateKMatrixMasstransport(void);
    … …  
    196197                ElementVector* CreatePVectorBalancethickness_CG(void);
    197198                ElementVector* CreatePVectorBalancevelocity(void);
     199                ElementVector* CreatePVectorSmoothedSlopeX(void);
     200                ElementVector* CreatePVectorSmoothedSlopeY(void);
    198201                ElementVector* CreatePVectorMasstransport(void);
    199202                ElementVector* CreatePVectorMasstransport_CG(void);

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r16007 r16010  
    301301        SurfaceSlopeXAnalysisEnum,
    302302        SurfaceSlopeYAnalysisEnum,
     303        SmoothedSurfaceSlopeXAnalysisEnum,
     304        SmoothedSurfaceSlopeYAnalysisEnum,
    303305        ThermalAnalysisEnum,
    304306        ThermalSolutionEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r16007 r16010  
    307307                case SurfaceSlopeXAnalysisEnum : return "SurfaceSlopeXAnalysis";
    308308                case SurfaceSlopeYAnalysisEnum : return "SurfaceSlopeYAnalysis";
     309                case SmoothedSurfaceSlopeXAnalysisEnum : return "SmoothedSurfaceSlopeXAnalysis";
     310                case SmoothedSurfaceSlopeYAnalysisEnum : return "SmoothedSurfaceSlopeYAnalysis";
    309311                case ThermalAnalysisEnum : return "ThermalAnalysis";
    310312                case ThermalSolutionEnum : return "ThermalSolution";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r16007 r16010  
    313313              else if (strcmp(name,"SurfaceSlopeXAnalysis")==0) return SurfaceSlopeXAnalysisEnum;
    314314              else if (strcmp(name,"SurfaceSlopeYAnalysis")==0) return SurfaceSlopeYAnalysisEnum;
     315              else if (strcmp(name,"SmoothedSurfaceSlopeXAnalysis")==0) return SmoothedSurfaceSlopeXAnalysisEnum;
     316              else if (strcmp(name,"SmoothedSurfaceSlopeYAnalysis")==0) return SmoothedSurfaceSlopeYAnalysisEnum;
    315317              else if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum;
    316318              else if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum;
    … …  
    381383              else if (strcmp(name,"Profiler")==0) return ProfilerEnum;
    382384              else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
    383               else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
    384               else if (strcmp(name,"VectorParam")==0) return VectorParamEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
     388              if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
     389              else if (strcmp(name,"VectorParam")==0) return VectorParamEnum;
     390              else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
    389391              else if (strcmp(name,"RiftfrontType")==0) return RiftfrontTypeEnum;
    390392              else if (strcmp(name,"Segment")==0) return SegmentEnum;
    … …  
    504506              else if (strcmp(name,"P1bubble")==0) return P1bubbleEnum;
    505507              else if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum;
    506               else if (strcmp(name,"P2")==0) return P2Enum;
    507               else if (strcmp(name,"P2xP1")==0) return P2xP1Enum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"P1xP2")==0) return P1xP2Enum;
     511              if (strcmp(name,"P2")==0) return P2Enum;
     512              else if (strcmp(name,"P2xP1")==0) return P2xP1Enum;
     513              else if (strcmp(name,"P1xP2")==0) return P1xP2Enum;
    512514              else if (strcmp(name,"P1P1")==0) return P1P1Enum;
    513515              else if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
Note: See TracChangeset for help on using the changeset viewer.