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.