Changeset 15877


Ignore:
Timestamp:
08/22/13 11:44:02 (12 years ago)
Author:
Mathieu Morlighem
Message:

NEW: added free surface solver for bottom and upper surface, not working yet

Location:
issm/trunk-jpl/src/c
Files:
10 added
16 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk-jpl/src/c/Makefile.am

    r15848 r15877  
    373373#Masstransport sources  {{{
    374374masstransport_sources = ./modules/ModelProcessorx/Masstransport/UpdateElementsMasstransport.cpp\
    375                                               ./modules/ModelProcessorx/Masstransport/CreateNodesMasstransport.cpp\
    376                                               ./modules/ModelProcessorx/Masstransport/CreateConstraintsMasstransport.cpp\
    377                                               ./modules/ModelProcessorx/Masstransport/CreateLoadsMasstransport.cpp\
    378                                                         ./analyses/masstransport_core.cpp
     375                                                                ./modules/ModelProcessorx/Masstransport/CreateNodesMasstransport.cpp\
     376                                                                ./modules/ModelProcessorx/Masstransport/CreateConstraintsMasstransport.cpp\
     377                                                                ./modules/ModelProcessorx/Masstransport/CreateLoadsMasstransport.cpp\
     378                                                                ./modules/ModelProcessorx/FreeSurfaceTop/UpdateElementsFreeSurfaceTop.cpp\
     379                                                                ./modules/ModelProcessorx/FreeSurfaceTop/CreateNodesFreeSurfaceTop.cpp\
     380                                                                ./modules/ModelProcessorx/FreeSurfaceTop/CreateConstraintsFreeSurfaceTop.cpp\
     381                                                                ./modules/ModelProcessorx/FreeSurfaceTop/CreateLoadsFreeSurfaceTop.cpp\
     382                                                                ./modules/ModelProcessorx/FreeSurfaceBase/UpdateElementsFreeSurfaceBase.cpp\
     383                                                                ./modules/ModelProcessorx/FreeSurfaceBase/CreateNodesFreeSurfaceBase.cpp\
     384                                                                ./modules/ModelProcessorx/FreeSurfaceBase/CreateConstraintsFreeSurfaceBase.cpp\
     385                                                                ./modules/ModelProcessorx/FreeSurfaceBase/CreateLoadsFreeSurfaceBase.cpp\
     386                                                                ./analyses/masstransport_core.cpp
     387
    379388#}}}
    380389#Thermal sources  {{{
  • issm/trunk-jpl/src/c/analyses/AnalysisConfiguration.cpp

    r15771 r15877  
    108108
    109109                case TransientSolutionEnum:
    110                         numanalyses=10-1;
     110                        numanalyses=11;
    111111                        analyses=xNew<int>(numanalyses);
    112                         analyses[0]=StressbalanceAnalysisEnum;
    113                         analyses[1]=StressbalanceVerticalAnalysisEnum;
    114                         analyses[2]=StressbalanceSIAAnalysisEnum;
    115                         analyses[3]=SurfaceSlopeAnalysisEnum;
    116                         analyses[4]=BedSlopeAnalysisEnum;
    117                         analyses[5]=ThermalAnalysisEnum;
    118                         analyses[6]=MeltingAnalysisEnum;
    119                         analyses[7]=EnthalpyAnalysisEnum;
    120                         analyses[8]=MasstransportAnalysisEnum;
     112                        analyses[ 0]=StressbalanceAnalysisEnum;
     113                        analyses[ 1]=StressbalanceVerticalAnalysisEnum;
     114                        analyses[ 2]=StressbalanceSIAAnalysisEnum;
     115                        analyses[ 3]=SurfaceSlopeAnalysisEnum;
     116                        analyses[ 4]=BedSlopeAnalysisEnum;
     117                        analyses[ 5]=ThermalAnalysisEnum;
     118                        analyses[ 6]=MeltingAnalysisEnum;
     119                        analyses[ 7]=EnthalpyAnalysisEnum;
     120                        analyses[ 8]=MasstransportAnalysisEnum;
     121                        analyses[ 9]=FreeSurfaceBaseAnalysisEnum;
     122                        analyses[10]=FreeSurfaceTopAnalysisEnum;
    121123                        break;
    122124
  • issm/trunk-jpl/src/c/analyses/masstransport_core.cpp

    r15849 r15877  
    1313
    1414        /*parameters: */
    15         bool save_results;
    16         bool issmbgradients,ispdd,isdelta18o;
    17         int  solution_type;
     15        bool  save_results;
     16        bool  issmbgradients,ispdd,isdelta18o,isFS,isfreesurface;
     17        int   solution_type;
    1818        int  *requested_outputs = NULL;
    19         int  numoutputs=0;
     19        int   numoutputs        = 0;
    2020
    21         /*activate formulation: */
     21        /*activate configuration*/
    2222        femmodel->SetCurrentConfiguration(MasstransportAnalysisEnum);
    2323
     
    2727        femmodel->parameters->FindParam(&ispdd,SurfaceforcingsIspddEnum);
    2828        femmodel->parameters->FindParam(&isdelta18o,SurfaceforcingsIsdelta18oEnum);
     29        femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
     30        femmodel->parameters->FindParam(&isfreesurface,MasstransportIsfreesurfaceEnum);
    2931        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    3032        femmodel->parameters->FindParam(&numoutputs,MasstransportNumRequestedOutputsEnum);
     
    4345                PositiveDegreeDayx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    4446        }
    45         if(VerboseSolution()) _printf0_("   call computational core\n");
    46         solutionsequence_linear(femmodel);
     47
     48        if(isFS && isfreesurface){
     49                if(VerboseSolution()) _printf0_("   call free surface computational core\n");
     50                femmodel->SetCurrentConfiguration(FreeSurfaceBaseAnalysisEnum);
     51                solutionsequence_linear(femmodel);
     52                femmodel->SetCurrentConfiguration(FreeSurfaceTopAnalysisEnum);
     53                solutionsequence_linear(femmodel);
     54        }
     55        else{
     56                if(VerboseSolution()) _printf0_("   call computational core\n");
     57                solutionsequence_linear(femmodel);
     58        }
    4759
    4860        if(save_results){
    4961                if(VerboseSolution()) _printf0_("   saving results\n");
    5062                InputToResultx(femmodel,ThicknessEnum);
     63                InputToResultx(femmodel,BedEnum);
     64                InputToResultx(femmodel,SurfaceEnum);
    5165                femmodel->RequestedOutputsx(requested_outputs,numoutputs);
    5266        }
  • issm/trunk-jpl/src/c/analyses/transient_core.cpp

    r15849 r15877  
    154154                if(save_results){
    155155                        if(VerboseSolution()) _printf0_("   saving transient results\n");
    156                         InputToResultx(femmodel,SurfaceEnum);
    157                         InputToResultx(femmodel,BedEnum);
    158156                        InputToResultx(femmodel,SurfaceforcingsMassBalanceEnum);
    159157                        InputToResultx(femmodel,MaskElementonfloatingiceEnum);
    160158                        femmodel->RequestedOutputsx(requested_outputs,numoutputs);
    161 
    162159                        if(isdelta18o){
    163160                                InputToResultx(femmodel,SurfaceforcingsMonthlytemperaturesEnum);
    164161                                InputToResultx(femmodel,SurfaceforcingsPrecipitationEnum);
    165                                 InputToResultx(femmodel,BasalFrictionEnum);
    166162                        }
    167163                        if(isgroundingline && (groundingline_migration==SubelementMigrationEnum || groundingline_migration==SubelementMigration2Enum)){
    168164                                InputToResultx(femmodel,GLlevelsetEnum);
    169165                        }
    170 
    171166                        if(VerboseSolution()) _printf0_("   saving temporary results\n");
    172167                        OutputResultsx(femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,femmodel->results);
  • issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

    r15872 r15877  
    469469                case MasstransportAnalysisEnum:
    470470                        return CreateKMatrixMasstransport();
     471                        break;
     472                case FreeSurfaceTopAnalysisEnum:
     473                        return CreateKMatrixFreeSurfaceTop();
     474                        break;
     475                case FreeSurfaceBaseAnalysisEnum:
     476                        return CreateKMatrixFreeSurfaceBase();
    471477                        break;
    472478                #ifdef _HAVE_BALANCED_
     
    518524        this->inputs->DeleteInput(VxAverageEnum);
    519525        this->inputs->DeleteInput(VyAverageEnum);
     526
     527        /*clean up and return*/
     528        return Ke;
     529}
     530/*}}}*/
     531/*FUNCTION Penta::CreateKMatrixFreeSurfaceTop {{{*/
     532ElementMatrix* Penta::CreateKMatrixFreeSurfaceTop(void){
     533
     534        if(!IsOnSurface()) return NULL;
     535
     536        Tria* tria=(Tria*)SpawnTria(1); //lower face is 0, upper face is 1.
     537        ElementMatrix* Ke=tria->CreateKMatrixFreeSurfaceTop();
     538        delete tria->material; delete tria;
     539
     540        /*clean up and return*/
     541        return Ke;
     542}
     543/*}}}*/
     544/*FUNCTION Penta::CreateKMatrixFreeSurfaceBase {{{*/
     545ElementMatrix* Penta::CreateKMatrixFreeSurfaceBase(void){
     546
     547        if(!IsOnBed()) return NULL;
     548
     549        Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
     550        ElementMatrix* Ke=tria->CreateKMatrixFreeSurfaceBase();
     551        delete tria->material; delete tria;
    520552
    521553        /*clean up and return*/
     
    666698                        return CreatePVectorMasstransport();
    667699                        break;
     700                case FreeSurfaceTopAnalysisEnum:
     701                        return CreatePVectorFreeSurfaceTop();
     702                        break;
     703                case FreeSurfaceBaseAnalysisEnum:
     704                        return CreatePVectorFreeSurfaceBase();
     705                        break;
    668706                #ifdef _HAVE_BALANCED_
    669707                case BalancethicknessAnalysisEnum:
     
    704742        this->inputs->DeleteInput(VxAverageEnum);
    705743        this->inputs->DeleteInput(VyAverageEnum);
     744
     745        /*Clean up and return*/
     746        return pe;
     747}
     748/*}}}*/
     749/*FUNCTION Penta::CreatePVectorFreeSurfaceTop {{{*/
     750ElementVector* Penta::CreatePVectorFreeSurfaceTop(void){
     751
     752        if(!IsOnSurface()) return NULL;
     753
     754        /*Call Tria function*/
     755        Tria* tria=(Tria*)SpawnTria(1); //lower face is 0, upper face is 1.
     756        ElementVector* pe=tria->CreatePVectorFreeSurfaceTop();
     757        delete tria->material; delete tria;
     758
     759        /*Clean up and return*/
     760        return pe;
     761}
     762/*}}}*/
     763/*FUNCTION Penta::CreatePVectorFreeSurfaceBase {{{*/
     764ElementVector* Penta::CreatePVectorFreeSurfaceBase(void){
     765
     766        if(!IsOnBed()) return NULL;
     767
     768        /*Call Tria function*/
     769        Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
     770        ElementVector* pe=tria->CreatePVectorFreeSurfaceBase();
     771        delete tria->material; delete tria;
    706772
    707773        /*Clean up and return*/
     
    22632329                InputUpdateFromSolutionMasstransport(solution);
    22642330                break;
     2331        case FreeSurfaceTopAnalysisEnum:
     2332                InputUpdateFromSolutionFreeSurfaceTop(solution);
     2333                break;
     2334        case FreeSurfaceBaseAnalysisEnum:
     2335                InputUpdateFromSolutionFreeSurfaceBase(solution);
     2336                break;
    22652337        #ifdef _HAVE_BALANCED_
    22662338        case BalancethicknessAnalysisEnum:
     
    23552427                /*Stop if we have reached the surface*/
    23562428                if (penta->IsOnSurface()) break;
     2429
     2430                /* get upper Penta*/
     2431                penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id);
     2432        }
     2433
     2434        /*Free ressources:*/
     2435        xDelete<int>(doflist);
     2436}
     2437/*}}}*/
     2438/*FUNCTION Penta::InputUpdateFromSolutionFreeSurfaceTop{{{*/
     2439void  Penta::InputUpdateFromSolutionFreeSurfaceTop(IssmDouble* solution){
     2440
     2441        const int  numdof   = NDOF1*NUMVERTICES;
     2442        const int  numdof2d = NDOF1*NUMVERTICES2D;
     2443
     2444        int    i;
     2445        int*   doflist = NULL;
     2446        IssmDouble newthickness[numdof];
     2447        IssmDouble newbed[numdof];
     2448        IssmDouble newsurface[numdof];
     2449
     2450        /*If not on bed, return*/
     2451        if (!IsOnSurface()) return;
     2452
     2453        /*Get dof list: */
     2454        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     2455
     2456        /*Use the dof list to index into the solution vector and extrude it */
     2457        for(i=0;i<numdof2d;i++){
     2458                newsurface[i+numdof2d]=solution[doflist[i+numdof2d]];
     2459                if(xIsNan<IssmDouble>(newsurface[i+numdof2d])) _error_("NaN found in solution vector");
     2460                newsurface[i]=newsurface[i+numdof2d];
     2461        }
     2462
     2463        /*Get previous bed and thickness*/
     2464        GetInputListOnVertices(&newbed[0],BedEnum);
     2465
     2466        for(i=0;i<numdof;i++) {
     2467                newthickness[i]=newsurface[i]-newbed[i];
     2468                _assert_(newthickness[i]>0.);
     2469        }
     2470
     2471        /*Start looping over all elements above current element and update all inputs*/
     2472        Penta* penta=this;
     2473        for(;;){
     2474                /*Add input to the element: */
     2475                penta->inputs->AddInput(new PentaInput(ThicknessEnum,newthickness,P1Enum));
     2476                penta->inputs->AddInput(new PentaInput(SurfaceEnum,newsurface,P1Enum));
     2477
     2478                /*Stop if we have reached the surface*/
     2479                if(penta->IsOnBed()) break;
     2480
     2481                /* get upper Penta*/
     2482                penta=penta->GetLowerElement(); _assert_(penta->Id()!=this->id);
     2483        }
     2484
     2485        /*Free ressources:*/
     2486        xDelete<int>(doflist);
     2487}
     2488/*}}}*/
     2489/*FUNCTION Penta::InputUpdateFromSolutionFreeSurfaceBase{{{*/
     2490void  Penta::InputUpdateFromSolutionFreeSurfaceBase(IssmDouble* solution){
     2491
     2492        const int  numdof   = NDOF1*NUMVERTICES;
     2493        const int  numdof2d = NDOF1*NUMVERTICES2D;
     2494
     2495        int    i;
     2496        int*   doflist = NULL;
     2497        IssmDouble newthickness[numdof];
     2498        IssmDouble newbed[numdof];
     2499        IssmDouble newsurface[numdof];
     2500
     2501        /*If not on bed, return*/
     2502        if (!IsOnBed()) return;
     2503
     2504        /*Get dof list: */
     2505        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     2506
     2507        /*Use the dof list to index into the solution vector and extrude it */
     2508        for(i=0;i<numdof2d;i++){
     2509                newbed[i]=solution[doflist[i]];
     2510                if(xIsNan<IssmDouble>(newbed[i])) _error_("NaN found in solution vector");
     2511                newbed[i+numdof2d]=newbed[i];
     2512        }
     2513
     2514        /*Get previous bed and thickness*/
     2515        GetInputListOnVertices(&newsurface[0],SurfaceEnum);
     2516
     2517        for(i=0;i<numdof;i++) {
     2518                newthickness[i]=newsurface[i]-newbed[i];
     2519                _assert_(newthickness[i]>0.);
     2520        }
     2521
     2522        /*Start looping over all elements above current element and update all inputs*/
     2523        Penta* penta=this;
     2524        for(;;){
     2525                /*Add input to the element: */
     2526                penta->inputs->AddInput(new PentaInput(ThicknessEnum,newthickness,P1Enum));
     2527                penta->inputs->AddInput(new PentaInput(BedEnum,newbed,P1Enum));
     2528
     2529                /*Stop if we have reached the surface*/
     2530                if(penta->IsOnSurface()) break;
    23572531
    23582532                /* get upper Penta*/
  • issm/trunk-jpl/src/c/classes/Elements/Penta.h

    r15867 r15877  
    181181                ElementMatrix* CreateKMatrix(void);
    182182                ElementMatrix* CreateKMatrixMasstransport(void);
     183                ElementMatrix* CreateKMatrixFreeSurfaceTop(void);
     184                ElementMatrix* CreateKMatrixFreeSurfaceBase(void);
    183185                ElementVector* CreatePVector(void);
    184186                ElementVector* CreatePVectorMasstransport(void);
     187                ElementVector* CreatePVectorFreeSurfaceTop(void);
     188                ElementVector* CreatePVectorFreeSurfaceBase(void);
    185189                ElementVector* CreatePVectorSlope(void);
    186190                void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[6][3],int numpoints);
     
    211215                void             InputExtrude(int enum_type,int object_type);
    212216                void           InputUpdateFromSolutionMasstransport(IssmDouble* solutiong);
     217                void           InputUpdateFromSolutionFreeSurfaceTop(IssmDouble* solutiong);
     218                void           InputUpdateFromSolutionFreeSurfaceBase(IssmDouble* solutiong);
    213219                void           InputUpdateFromSolutionOneDof(IssmDouble* solutiong,int enum_type);
    214220                void           InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solutiong,int enum_type);
  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r15814 r15877  
    64876487}
    64886488/*}}}*/
     6489/*FUNCTION Tria::CreateKMatrixFreeSurfaceTop {{{*/
     6490ElementMatrix* Tria::CreateKMatrixFreeSurfaceTop(void){
     6491
     6492        /*Intermediaries */
     6493        int        stabilization;
     6494        int        dim;
     6495        IssmDouble Jdettria,D_scalar,dt,h;
     6496        IssmDouble vel,vx,vy,dvxdx,dvydy;
     6497        IssmDouble dvx[2],dvy[2];
     6498        IssmDouble xyz_list[NUMVERTICES][3];
     6499
     6500        /*Fetch number of nodes for this finite element*/
     6501        int numnodes = this->NumberofNodes();
     6502
     6503        /*Initialize Element matrix and vectors*/
     6504        ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
     6505        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     6506        IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
     6507        IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
     6508        IssmDouble     D[2][2];
     6509
     6510        /*Retrieve all inputs and parameters*/
     6511        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     6512        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     6513        this->parameters->FindParam(&stabilization,MasstransportStabilizationEnum);
     6514        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     6515        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
     6516        h=sqrt(2*this->GetArea());
     6517
     6518        /* Start  looping on the number of gaussian points: */
     6519        GaussTria *gauss=new GaussTria(2);
     6520        for(int ig=gauss->begin();ig<gauss->end();ig++){
     6521
     6522                gauss->GaussPoint(ig);
     6523
     6524                GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
     6525                GetNodalFunctions(basis,gauss);
     6526
     6527                vx_input->GetInputValue(&vx,gauss);
     6528                vy_input->GetInputValue(&vy,gauss);
     6529                vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
     6530                vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
     6531
     6532                D_scalar=gauss->weight*Jdettria;
     6533
     6534                TripleMultiply(basis,1,numnodes,1,
     6535                                        &D_scalar,1,1,0,
     6536                                        basis,1,numnodes,0,
     6537                                        &Ke->values[0],1);
     6538
     6539                GetBMasstransport(B,&xyz_list[0][0],gauss);
     6540                GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);
     6541
     6542                dvxdx=dvx[0];
     6543                dvydy=dvy[1];
     6544                D_scalar=dt*gauss->weight*Jdettria;
     6545
     6546                D[0][0]=D_scalar*vx;
     6547                D[0][1]=0.;
     6548                D[1][0]=0.;
     6549                D[1][1]=D_scalar*vy;
     6550                TripleMultiply(B,2,numnodes,1,
     6551                                        &D[0][0],2,2,0,
     6552                                        Bprime,2,numnodes,0,
     6553                                        &Ke->values[0],1);
     6554
     6555                if(stabilization==2){
     6556                        /*Streamline upwinding*/
     6557                        vel=sqrt(vx*vx+vy*vy)+1.e-8;
     6558                        D[0][0]=h/(2*vel)*vx*vx;
     6559                        D[1][0]=h/(2*vel)*vy*vx;
     6560                        D[0][1]=h/(2*vel)*vx*vy;
     6561                        D[1][1]=h/(2*vel)*vy*vy;
     6562                }
     6563                else if(stabilization==1){
     6564                        /*SSA*/
     6565                        vx_input->GetInputAverage(&vx);
     6566                        vy_input->GetInputAverage(&vy);
     6567                        D[0][0]=h/2.0*fabs(vx);
     6568                        D[0][1]=0.;
     6569                        D[1][0]=0.;
     6570                        D[1][1]=h/2.0*fabs(vy);
     6571                }
     6572                if(stabilization==1 || stabilization==2){
     6573                        D[0][0]=D_scalar*D[0][0];
     6574                        D[1][0]=D_scalar*D[1][0];
     6575                        D[0][1]=D_scalar*D[0][1];
     6576                        D[1][1]=D_scalar*D[1][1];
     6577                        TripleMultiply(Bprime,2,numnodes,1,
     6578                                                &D[0][0],2,2,0,
     6579                                                Bprime,2,numnodes,0,
     6580                                                &Ke->values[0],1);
     6581                }
     6582        }
     6583
     6584        /*Clean up and return*/
     6585        xDelete<IssmDouble>(basis);
     6586        xDelete<IssmDouble>(B);
     6587        xDelete<IssmDouble>(Bprime);
     6588        delete gauss;
     6589        return Ke;
     6590}
     6591/*}}}*/
     6592/*FUNCTION Tria::CreateKMatrixFreeSurfaceBase {{{*/
     6593ElementMatrix* Tria::CreateKMatrixFreeSurfaceBase(void){
     6594
     6595        /*Intermediaries */
     6596        int        stabilization;
     6597        int        dim;
     6598        IssmDouble Jdettria,D_scalar,dt,h;
     6599        IssmDouble vel,vx,vy,dvxdx,dvydy;
     6600        IssmDouble dvx[2],dvy[2];
     6601        IssmDouble xyz_list[NUMVERTICES][3];
     6602
     6603        /*Fetch number of nodes for this finite element*/
     6604        int numnodes = this->NumberofNodes();
     6605
     6606        /*Initialize Element matrix and vectors*/
     6607        ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
     6608        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     6609        IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
     6610        IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
     6611        IssmDouble     D[2][2];
     6612
     6613        /*Retrieve all inputs and parameters*/
     6614        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     6615        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     6616        this->parameters->FindParam(&stabilization,MasstransportStabilizationEnum);
     6617        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     6618        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
     6619        h=sqrt(2*this->GetArea());
     6620
     6621        /* Start  looping on the number of gaussian points: */
     6622        GaussTria *gauss=new GaussTria(2);
     6623        for(int ig=gauss->begin();ig<gauss->end();ig++){
     6624
     6625                gauss->GaussPoint(ig);
     6626
     6627                GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
     6628                GetNodalFunctions(basis,gauss);
     6629
     6630                vx_input->GetInputValue(&vx,gauss);
     6631                vy_input->GetInputValue(&vy,gauss);
     6632                vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
     6633                vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
     6634
     6635                D_scalar=gauss->weight*Jdettria;
     6636
     6637                TripleMultiply(basis,1,numnodes,1,
     6638                                        &D_scalar,1,1,0,
     6639                                        basis,1,numnodes,0,
     6640                                        &Ke->values[0],1);
     6641
     6642                GetBMasstransport(B,&xyz_list[0][0],gauss);
     6643                GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);
     6644
     6645                dvxdx=dvx[0];
     6646                dvydy=dvy[1];
     6647                D_scalar=dt*gauss->weight*Jdettria;
     6648
     6649                D[0][0]=D_scalar*vx;
     6650                D[0][1]=0.;
     6651                D[1][0]=0.;
     6652                D[1][1]=D_scalar*vy;
     6653                TripleMultiply(B,2,numnodes,1,
     6654                                        &D[0][0],2,2,0,
     6655                                        Bprime,2,numnodes,0,
     6656                                        &Ke->values[0],1);
     6657
     6658                if(stabilization==2){
     6659                        /*Streamline upwinding*/
     6660                        vel=sqrt(vx*vx+vy*vy)+1.e-8;
     6661                        D[0][0]=h/(2*vel)*vx*vx;
     6662                        D[1][0]=h/(2*vel)*vy*vx;
     6663                        D[0][1]=h/(2*vel)*vx*vy;
     6664                        D[1][1]=h/(2*vel)*vy*vy;
     6665                }
     6666                else if(stabilization==1){
     6667                        /*SSA*/
     6668                        vx_input->GetInputAverage(&vx);
     6669                        vy_input->GetInputAverage(&vy);
     6670                        D[0][0]=h/2.0*fabs(vx);
     6671                        D[0][1]=0.;
     6672                        D[1][0]=0.;
     6673                        D[1][1]=h/2.0*fabs(vy);
     6674                }
     6675                if(stabilization==1 || stabilization==2){
     6676                        D[0][0]=D_scalar*D[0][0];
     6677                        D[1][0]=D_scalar*D[1][0];
     6678                        D[0][1]=D_scalar*D[0][1];
     6679                        D[1][1]=D_scalar*D[1][1];
     6680                        TripleMultiply(Bprime,2,numnodes,1,
     6681                                                &D[0][0],2,2,0,
     6682                                                Bprime,2,numnodes,0,
     6683                                                &Ke->values[0],1);
     6684                }
     6685        }
     6686
     6687        /*Clean up and return*/
     6688        xDelete<IssmDouble>(basis);
     6689        xDelete<IssmDouble>(B);
     6690        xDelete<IssmDouble>(Bprime);
     6691        delete gauss;
     6692        return Ke;
     6693}
     6694/*}}}*/
    64896695/*FUNCTION Tria::CreatePVectorMasstransport{{{*/
    64906696ElementVector* Tria::CreatePVectorMasstransport(void){
     
    65056711        /*Intermediaries */
    65066712        IssmDouble Jdettria,dt;
    6507         IssmDouble surface_mass_balance_g,basal_melting_g,basal_melting_correction_g,thickness_g;
     6713        IssmDouble ms,mb,mb_correction,thickness;
    65086714        IssmDouble xyz_list[NUMVERTICES][3];
    65096715
     
    65186724        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    65196725        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    6520         Input* surface_mass_balance_input     = inputs->GetInput(SurfaceforcingsMassBalanceEnum);         _assert_(surface_mass_balance_input);
    6521         Input* basal_melting_input            = inputs->GetInput(BasalforcingsMeltingRateEnum);           _assert_(basal_melting_input);
    6522         Input* basal_melting_correction_input = inputs->GetInput(BasalforcingsMeltingRateCorrectionEnum);
    6523         Input* thickness_input                = inputs->GetInput(ThicknessEnum);                          _assert_(thickness_input);
    6524 
    6525         /*Initialize basal_melting_correction_g to 0, do not forget!:*/
     6726        Input* ms_input     = inputs->GetInput(SurfaceforcingsMassBalanceEnum);      _assert_(ms_input);
     6727        Input* mb_input     = inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(mb_input);
     6728        Input* mb_correction_input = inputs->GetInput(BasalforcingsMeltingRateCorrectionEnum);
     6729        Input* thickness_input  = inputs->GetInput(ThicknessEnum);     _assert_(thickness_input);
     6730
     6731        /*Initialize mb_correction to 0, do not forget!:*/
    65266732        /* Start  looping on the number of gaussian points: */
    65276733        GaussTria* gauss=new GaussTria(2);
     
    65336739                GetNodalFunctions(basis,gauss);
    65346740
    6535                 surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
    6536                 basal_melting_input->GetInputValue(&basal_melting_g,gauss);
    6537                 thickness_input->GetInputValue(&thickness_g,gauss);
    6538                 if(basal_melting_correction_input)
    6539                  basal_melting_correction_input->GetInputValue(&basal_melting_correction_g,gauss);
     6741                ms_input->GetInputValue(&ms,gauss);
     6742                mb_input->GetInputValue(&mb,gauss);
     6743                thickness_input->GetInputValue(&thickness,gauss);
     6744                if(mb_correction_input)
     6745                 mb_correction_input->GetInputValue(&mb_correction,gauss);
    65406746                else
    6541                  basal_melting_correction_g=0.;
    6542 
    6543                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g-basal_melting_correction_g))*basis[i];
     6747                 mb_correction=0.;
     6748
     6749                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness+dt*(ms-mb-mb_correction))*basis[i];
    65446750        }
    65456751
     
    65556761        /*Intermediaries */
    65566762        IssmDouble Jdettria,dt;
    6557         IssmDouble surface_mass_balance_g,basal_melting_g,thickness_g;
     6763        IssmDouble ms,mb,thickness;
    65586764        IssmDouble xyz_list[NUMVERTICES][3];
    65596765
     
    65686774        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    65696775        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6570         Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input);
    6571         Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum);          _assert_(basal_melting_input);
    6572         Input* thickness_input=inputs->GetInput(ThicknessEnum);                             _assert_(thickness_input);
     6776        Input* ms_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(ms_input);
     6777        Input* mb_input=inputs->GetInput(BasalforcingsMeltingRateEnum);   _assert_(mb_input);
     6778        Input* thickness_input=inputs->GetInput(ThicknessEnum);           _assert_(thickness_input);
    65736779
    65746780        /* Start  looping on the number of gaussian points: */
     
    65816787                GetNodalFunctions(basis,gauss);
    65826788
    6583                 surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
    6584                 basal_melting_input->GetInputValue(&basal_melting_g,gauss);
    6585                 thickness_input->GetInputValue(&thickness_g,gauss);
    6586 
    6587                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g))*basis[i];
     6789                ms_input->GetInputValue(&ms,gauss);
     6790                mb_input->GetInputValue(&mb,gauss);
     6791                thickness_input->GetInputValue(&thickness,gauss);
     6792
     6793                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness+dt*(ms-mb))*basis[i];
     6794        }
     6795
     6796        /*Clean up and return*/
     6797        xDelete<IssmDouble>(basis);
     6798        delete gauss;
     6799        return pe;
     6800}
     6801/*}}}*/
     6802/*FUNCTION Tria::CreatePVectorFreeSurfaceTop {{{*/
     6803ElementVector* Tria::CreatePVectorFreeSurfaceTop(void){
     6804
     6805        /*Intermediaries */
     6806        IssmDouble Jdettria,dt;
     6807        IssmDouble ms,surface,vz;
     6808        IssmDouble xyz_list[NUMVERTICES][3];
     6809
     6810        /*Fetch number of nodes and dof for this finite element*/
     6811        int numnodes = this->NumberofNodes();
     6812
     6813        /*Initialize Element vector and other vectors*/
     6814        ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
     6815        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     6816
     6817        /*Retrieve all inputs and parameters*/
     6818        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     6819        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     6820        Input* vz_input     = inputs->GetInput(VzEnum);                         _assert_(vz_input);
     6821        Input* ms_input     = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(ms_input);
     6822        Input* surface_input= inputs->GetInput(ThicknessEnum);              _assert_(surface_input);
     6823
     6824        /*Initialize mb_correction to 0, do not forget!:*/
     6825        /* Start  looping on the number of gaussian points: */
     6826        GaussTria* gauss=new GaussTria(2);
     6827        for(int ig=gauss->begin();ig<gauss->end();ig++){
     6828
     6829                gauss->GaussPoint(ig);
     6830
     6831                GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
     6832                GetNodalFunctions(basis,gauss);
     6833
     6834                ms_input->GetInputValue(&ms,gauss);
     6835                surface_input->GetInputValue(&surface,gauss);
     6836
     6837                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(surface + dt*ms + dt*vz)*basis[i];
     6838        }
     6839
     6840        /*Clean up and return*/
     6841        xDelete<IssmDouble>(basis);
     6842        delete gauss;
     6843        return pe;
     6844}
     6845/*}}}*/
     6846/*FUNCTION Tria::CreatePVectorFreeSurfaceBase {{{*/
     6847ElementVector* Tria::CreatePVectorFreeSurfaceBase(void){
     6848
     6849        /*Intermediaries */
     6850        IssmDouble Jdettria,dt;
     6851        IssmDouble mb,mb_correction,bed,vz;
     6852        IssmDouble xyz_list[NUMVERTICES][3];
     6853
     6854        /*Fetch number of nodes and dof for this finite element*/
     6855        int numnodes = this->NumberofNodes();
     6856
     6857        /*Initialize Element vector and other vectors*/
     6858        ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
     6859        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     6860
     6861        /*Retrieve all inputs and parameters*/
     6862        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     6863        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     6864        Input* vz_input  = inputs->GetInput(VzEnum);                         _assert_(vz_input);
     6865        Input* mb_input  = inputs->GetInput(BasalforcingsMeltingRateEnum);   _assert_(mb_input);
     6866        Input* mb_correction_input = inputs->GetInput(BasalforcingsMeltingRateCorrectionEnum);
     6867        Input* bed_input = inputs->GetInput(BedEnum);                        _assert_(bed_input);
     6868
     6869        /*Initialize mb_correction to 0, do not forget!:*/
     6870        /* Start  looping on the number of gaussian points: */
     6871        GaussTria* gauss=new GaussTria(2);
     6872        for(int ig=gauss->begin();ig<gauss->end();ig++){
     6873
     6874                gauss->GaussPoint(ig);
     6875
     6876                GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
     6877                GetNodalFunctions(basis,gauss);
     6878
     6879                mb_input->GetInputValue(&mb,gauss);
     6880                bed_input->GetInputValue(&bed,gauss);
     6881                if(mb_correction_input)
     6882                 mb_correction_input->GetInputValue(&mb_correction,gauss);
     6883                else
     6884                 mb_correction=0.;
     6885
     6886                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(bed+dt*(mb-mb_correction) + dt*vz)*basis[i];
    65886887        }
    65896888
     
    69387237        /*Intermediaries */
    69397238        IssmDouble xyz_list[NUMVERTICES][3];
    6940         IssmDouble dhdt_g,basal_melting_g,surface_mass_balance_g,Jdettria;
     7239        IssmDouble dhdt_g,mb_g,ms_g,Jdettria;
    69417240
    69427241        /*Fetch number of nodes and dof for this finite element*/
     
    69497248        /*Retrieve all inputs and parameters*/
    69507249        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6951         Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input);
    6952         Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum);          _assert_(basal_melting_input);
    6953         Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum);             _assert_(dhdt_input);
     7250        Input* ms_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(ms_input);
     7251        Input* mb_input=inputs->GetInput(BasalforcingsMeltingRateEnum);   _assert_(mb_input);
     7252        Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input);
    69547253
    69557254        /* Start  looping on the number of gaussian points: */
     
    69597258                gauss->GaussPoint(ig);
    69607259
    6961                 surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
    6962                 basal_melting_input->GetInputValue(&basal_melting_g,gauss);
     7260                ms_input->GetInputValue(&ms_g,gauss);
     7261                mb_input->GetInputValue(&mb_g,gauss);
    69637262                dhdt_input->GetInputValue(&dhdt_g,gauss);
    69647263
     
    69667265                GetNodalFunctions(basis,gauss);
    69677266
    6968                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*basis[i];
     7267                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(ms_g-mb_g-dhdt_g)*basis[i];
    69697268        }
    69707269
     
    69807279        /*Intermediaries */
    69817280        IssmDouble xyz_list[NUMVERTICES][3];
    6982         IssmDouble basal_melting_g,surface_mass_balance_g,dhdt_g,Jdettria;
     7281        IssmDouble mb_g,ms_g,dhdt_g,Jdettria;
    69837282
    69847283        /*Fetch number of nodes and dof for this finite element*/
     
    69917290        /*Retrieve all inputs and parameters*/
    69927291        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6993         Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input);
    6994         Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum);          _assert_(basal_melting_input);
    6995         Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum);                                       _assert_(dhdt_input);
     7292        Input* ms_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(ms_input);
     7293        Input* mb_input=inputs->GetInput(BasalforcingsMeltingRateEnum);   _assert_(mb_input);
     7294        Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum);_assert_(dhdt_input);
    69967295
    69977296        /* Start  looping on the number of gaussian points: */
     
    70017300                gauss->GaussPoint(ig);
    70027301
    7003                 surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
    7004                 basal_melting_input->GetInputValue(&basal_melting_g,gauss);
     7302                ms_input->GetInputValue(&ms_g,gauss);
     7303                mb_input->GetInputValue(&mb_g,gauss);
    70057304                dhdt_input->GetInputValue(&dhdt_g,gauss);
    70067305
     
    70087307                GetNodalFunctions(basis,gauss);
    70097308
    7010                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*basis[i];
     7309                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(ms_g-mb_g-dhdt_g)*basis[i];
    70117310        }
    70127311
  • issm/trunk-jpl/src/c/classes/Elements/Tria.h

    r15866 r15877  
    188188                ElementMatrix* CreateKMatrixMasstransport_CG(void);
    189189                ElementMatrix* CreateKMatrixMasstransport_DG(void);
     190                ElementMatrix* CreateKMatrixFreeSurfaceTop(void);
     191                ElementMatrix* CreateKMatrixFreeSurfaceBase(void);
    190192                ElementMatrix* CreateMassMatrix(void);
    191193                ElementVector* CreatePVector(void);
     
    196198                ElementVector* CreatePVectorMasstransport_CG(void);
    197199                ElementVector* CreatePVectorMasstransport_DG(void);
     200                ElementVector* CreatePVectorFreeSurfaceTop(void);
     201                ElementVector* CreatePVectorFreeSurfaceBase(void);
    198202                ElementVector* CreatePVectorSlope(void);
    199203                IssmDouble     GetArea(void);
  • issm/trunk-jpl/src/c/classes/Node.cpp

    r15862 r15877  
    9393        /*2d solutions in 3d, we need to constrain all the nodes that are not on base*/
    9494        if(
     95                                analysis_type==FreeSurfaceBaseAnalysisEnum ||
    9596                                analysis_type==MasstransportAnalysisEnum ||
    9697                                analysis_type==MeltingAnalysisEnum ||
     
    105106                        _assert_(iomodel->Data(MeshVertexonbedEnum));
    106107                        if(!(reCast<bool>(iomodel->Data(MeshVertexonbedEnum)[io_index]))){
     108                                this->Deactivate();
     109                        }
     110                }
     111        }
     112        if(analysis_type==FreeSurfaceTopAnalysisEnum){
     113                if(iomodel->dim==3){
     114                        /*On a 3d mesh, we may have collapsed elements, hence dead nodes. Freeze them out: */
     115                        _assert_(iomodel->Data(MeshVertexonsurfaceEnum));
     116                        if(!(reCast<bool>(iomodel->Data(MeshVertexonsurfaceEnum)[io_index]))){
    107117                                this->Deactivate();
    108118                        }
  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateDataSets.cpp

    r15771 r15877  
    148148                        UpdateElementsMasstransport(elements,iomodel,analysis_counter,analysis_type);
    149149                        break;
     150                case FreeSurfaceTopAnalysisEnum:
     151                        CreateNodesFreeSurfaceTop(pnodes, iomodel);
     152                        CreateConstraintsFreeSurfaceTop(pconstraints,iomodel);
     153                        CreateLoadsFreeSurfaceTop(ploads,iomodel);
     154                        UpdateElementsFreeSurfaceTop(elements,iomodel,analysis_counter,analysis_type);
     155                        break;
     156                case FreeSurfaceBaseAnalysisEnum:
     157                        CreateNodesFreeSurfaceBase(pnodes, iomodel);
     158                        CreateConstraintsFreeSurfaceBase(pconstraints,iomodel);
     159                        CreateLoadsFreeSurfaceBase(ploads,iomodel);
     160                        UpdateElementsFreeSurfaceBase(elements,iomodel,analysis_counter,analysis_type);
     161                        break;
    150162                #endif
    151163
  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r15771 r15877  
    5757        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceIsnewtonEnum));
    5858        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceMaxiterEnum));
     59        parameters->AddObject(iomodel->CopyConstantObject(StressbalancePenaltyFactorEnum));
    5960        parameters->AddObject(iomodel->CopyConstantObject(SteadystateReltolEnum));
    6061        parameters->AddObject(iomodel->CopyConstantObject(SteadystateMaxiterEnum));
     
    6566        parameters->AddObject(iomodel->CopyConstantObject(TimesteppingTimeStepEnum));
    6667        parameters->AddObject(iomodel->CopyConstantObject(TimesteppingCflCoefficientEnum));
     68        parameters->AddObject(iomodel->CopyConstantObject(MasstransportIsfreesurfaceEnum));
    6769        parameters->AddObject(iomodel->CopyConstantObject(MasstransportHydrostaticAdjustmentEnum));
    6870        parameters->AddObject(iomodel->CopyConstantObject(MasstransportStabilizationEnum));
    69         parameters->AddObject(iomodel->CopyConstantObject(StressbalancePenaltyFactorEnum));
    7071        parameters->AddObject(iomodel->CopyConstantObject(MasstransportMinThicknessEnum));
    7172        parameters->AddObject(iomodel->CopyConstantObject(MasstransportPenaltyFactorEnum));
  • issm/trunk-jpl/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp

    r15771 r15877  
    103103                        numdofs=1;
    104104                        break;
     105                case FreeSurfaceTopAnalysisEnum:
     106                        numdofs=1;
     107                        break;
     108                case FreeSurfaceBaseAnalysisEnum:
     109                        numdofs=1;
     110                        break;
    105111                case GiaAnalysisEnum:
    106112                        numdofs=1;
  • issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h

    r15771 r15877  
    104104void CreateLoadsMasstransport(Loads** ploads, IoModel* iomodel);
    105105void UpdateElementsMasstransport(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     106void CreateNodesFreeSurfaceTop(Nodes** pnodes,IoModel* iomodel);
     107void CreateConstraintsFreeSurfaceTop(Constraints** pconstraints,IoModel* iomodel);
     108void CreateLoadsFreeSurfaceTop(Loads** ploads, IoModel* iomodel);
     109void UpdateElementsFreeSurfaceTop(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     110void CreateNodesFreeSurfaceBase(Nodes** pnodes,IoModel* iomodel);
     111void CreateConstraintsFreeSurfaceBase(Constraints** pconstraints,IoModel* iomodel);
     112void CreateLoadsFreeSurfaceBase(Loads** ploads, IoModel* iomodel);
     113void UpdateElementsFreeSurfaceBase(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
    106114
    107115/*balancedthickness:*/
  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r15827 r15877  
    188188        MiscellaneousNameEnum, //FIXME: only used by qmu, should not be marshalled (already in queueing script)
    189189        MasstransportHydrostaticAdjustmentEnum,
     190        MasstransportIsfreesurfaceEnum,
    190191        MasstransportMinThicknessEnum,
    191192        MasstransportPenaltyFactorEnum,
     
    295296        MasstransportAnalysisEnum,
    296297        MasstransportSolutionEnum,
     298        FreeSurfaceBaseAnalysisEnum,
     299        FreeSurfaceTopAnalysisEnum,
    297300        SteadystateSolutionEnum,
    298301        SurfaceSlopeAnalysisEnum,
  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r15771 r15877  
    196196                case MiscellaneousNameEnum : return "MiscellaneousName";
    197197                case MasstransportHydrostaticAdjustmentEnum : return "MasstransportHydrostaticAdjustment";
     198                case MasstransportIsfreesurfaceEnum : return "MasstransportIsfreesurface";
    198199                case MasstransportMinThicknessEnum : return "MasstransportMinThickness";
    199200                case MasstransportPenaltyFactorEnum : return "MasstransportPenaltyFactor";
     
    301302                case MasstransportAnalysisEnum : return "MasstransportAnalysis";
    302303                case MasstransportSolutionEnum : return "MasstransportSolution";
     304                case FreeSurfaceBaseAnalysisEnum : return "FreeSurfaceBaseAnalysis";
     305                case FreeSurfaceTopAnalysisEnum : return "FreeSurfaceTopAnalysis";
    303306                case SteadystateSolutionEnum : return "SteadystateSolution";
    304307                case SurfaceSlopeAnalysisEnum : return "SurfaceSlopeAnalysis";
     
    587590                case PatersonEnum : return "Paterson";
    588591                case ArrheniusEnum : return "Arrhenius";
     592                case LliboutryDuvalEnum : return "LliboutryDuval";
    589593                case MaximumNumberOfDefinitionsEnum : return "MaximumNumberOfDefinitions";
    590594                default : return "unknown";
  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r15771 r15877  
    199199              else if (strcmp(name,"MiscellaneousName")==0) return MiscellaneousNameEnum;
    200200              else if (strcmp(name,"MasstransportHydrostaticAdjustment")==0) return MasstransportHydrostaticAdjustmentEnum;
     201              else if (strcmp(name,"MasstransportIsfreesurface")==0) return MasstransportIsfreesurfaceEnum;
    201202              else if (strcmp(name,"MasstransportMinThickness")==0) return MasstransportMinThicknessEnum;
    202203              else if (strcmp(name,"MasstransportPenaltyFactor")==0) return MasstransportPenaltyFactorEnum;
     
    259260              else if (strcmp(name,"GiaMantleViscosity")==0) return GiaMantleViscosityEnum;
    260261              else if (strcmp(name,"GiaLithosphereThickness")==0) return GiaLithosphereThicknessEnum;
    261               else if (strcmp(name,"Thickness")==0) return ThicknessEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"TimesteppingStartTime")==0) return TimesteppingStartTimeEnum;
     265              if (strcmp(name,"Thickness")==0) return ThicknessEnum;
     266              else if (strcmp(name,"TimesteppingStartTime")==0) return TimesteppingStartTimeEnum;
    266267              else if (strcmp(name,"TimesteppingFinalTime")==0) return TimesteppingFinalTimeEnum;
    267268              else if (strcmp(name,"TimesteppingCflCoefficient")==0) return TimesteppingCflCoefficientEnum;
     
    307308              else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum;
    308309              else if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum;
     310              else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum;
     311              else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum;
    309312              else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum;
    310313              else if (strcmp(name,"SurfaceSlopeAnalysis")==0) return SurfaceSlopeAnalysisEnum;
     
    380383              else if (strcmp(name,"Profiler")==0) return ProfilerEnum;
    381384              else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
    382               else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
    383               else if (strcmp(name,"VectorParam")==0) return VectorParamEnum;
    384               else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"RiftfrontType")==0) return RiftfrontTypeEnum;
     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;
     391              else if (strcmp(name,"RiftfrontType")==0) return RiftfrontTypeEnum;
    389392              else if (strcmp(name,"Segment")==0) return SegmentEnum;
    390393              else if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum;
     
    503506              else if (strcmp(name,"P1bubble")==0) return P1bubbleEnum;
    504507              else if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum;
    505               else if (strcmp(name,"P2")==0) return P2Enum;
    506               else if (strcmp(name,"P2xP1")==0) return P2xP1Enum;
    507               else if (strcmp(name,"P1xP2")==0) return P1xP2Enum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"P1P1")==0) return P1P1Enum;
     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;
     514              else if (strcmp(name,"P1P1")==0) return P1P1Enum;
    512515              else if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
    513516              else if (strcmp(name,"MINI")==0) return MINIEnum;
     
    599602              else if (strcmp(name,"Paterson")==0) return PatersonEnum;
    600603              else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum;
     604              else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum;
    601605              else if (strcmp(name,"MaximumNumberOfDefinitions")==0) return MaximumNumberOfDefinitionsEnum;
    602606         else stage=6;
Note: See TracChangeset for help on using the changeset viewer.