Changeset 8410


Ignore:
Timestamp:
05/24/11 11:29:04 (14 years ago)
Author:
seroussi
Message:

moved thermal PVectors from Tria to Penta

Location:
issm/trunk/src/c
Files:
6 edited
1 moved

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/Makefile.am

    r8404 r8410  
    635635                                        ./modules/ConstraintsStatex/RiftConstraintsState.cpp\
    636636                                        ./modules/ConstraintsStatex/ThermalConstraintsState.cpp\
    637                                         ./modules/ConstraintsStatex/MeltingIsPresent.cpp\
     637                                        ./modules/ConstraintsStatex/ThermalIsPresent.cpp\
    638638                                        ./modules/Responsex/Responsex.h\
    639639                                        ./modules/Responsex/Responsex.cpp\
    … …  
    12701270                                        ./modules/ConstraintsStatex/RiftConstraintsState.cpp\
    12711271                                        ./modules/ConstraintsStatex/ThermalConstraintsState.cpp\
    1272                                         ./modules/ConstraintsStatex/MeltingIsPresent.cpp\
     1272                                        ./modules/ConstraintsStatex/ThermalIsPresent.cpp\
    12731273                                        ./modules/Responsex/Responsex.h\
    12741274                                        ./modules/Responsex/Responsex.cpp\

  • issm/trunk/src/c/modules/ConstraintsStatex/ConstraintsStateLocal.h

    r8404 r8410  
    1111/*melting: */
    1212void  ThermalConstraintsState(Loads* loads, int* pconverged, int* pnum_unstable_constraints,int analysis_type);
    13 int   MeltingIsPresent(Loads* loads,int analysis_type);
     13int   ThermalIsPresent(Loads* loads,int analysis_type);
    1414
    1515/*rifts module: */

  • issm/trunk/src/c/modules/ConstraintsStatex/ConstraintsStatex.cpp

    r8404 r8410  
    3535                RiftConstraintsState(&converged,&num_unstable_constraints,loads,min_mechanical_constraints,analysis_type);
    3636        }
    37         else if(MeltingIsPresent(loads,analysis_type)){
     37        else if(ThermalIsPresent(loads,analysis_type)){
    3838                ThermalConstraintsState(loads,&converged,&num_unstable_constraints,analysis_type);
    3939        }

  • issm/trunk/src/c/modules/ConstraintsStatex/ThermalIsPresent.cpp

    r8407 r8410  
    1 /*!\file:  MeltingIsPresent.cpp
     1/*!\file:  ThermalIsPresent.cpp
    22 * \brief
    33 */
    … …  
    1111#include "./ConstraintsStateLocal.h"
    1212
    13 int MeltingIsPresent(Loads* loads,int configuration_type){
     13int ThermalIsPresent(Loads* loads,int configuration_type){
    1414
    1515        int i;

  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r8409 r8410  
    29842984ElementVector* Penta::CreatePVectorThermalShelf(void){
    29852985
     2986        /*Constants*/
     2987        const int  numdof=NUMVERTICES*NDOF1;
     2988
     2989        /*Intermediaries */
     2990        int        i,j,ig;
     2991        double     Jdet2d;
     2992        double     mixed_layer_capacity,thermal_exchange_velocity;
     2993        double     rho_ice,rho_water,pressure,dt,scalar_ocean;
     2994        double     meltingpoint,beta,heatcapacity,t_pmp;
     2995        double     xyz_list[NUMVERTICES][3];
     2996        double     xyz_list_tria[NUMVERTICES2D][3];
     2997        double     l1l6[NUMVERTICES];
     2998        GaussPenta* gauss=NULL;
     2999
     3000        /*Initialize Element vector*/
     3001        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     3002
    29863003        /* Ice/ocean heat exchange flux on ice shelf base */
    29873004        if (!IsOnBed() || !IsOnShelf()) return NULL;
    29883005
    2989         /*Call Tria function*/
    2990         Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
    2991         ElementVector* pe=tria->CreatePVectorThermalShelf();
    2992         delete tria->matice; delete tria;
     3006        /*Retrieve all inputs and parameters*/
     3007        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     3008        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
     3009        mixed_layer_capacity=matpar->GetMixedLayerCapacity();
     3010        thermal_exchange_velocity=matpar->GetThermalExchangeVelocity();
     3011        rho_water=matpar->GetRhoWater();
     3012        rho_ice=matpar->GetRhoIce();
     3013        heatcapacity=matpar->GetHeatCapacity();
     3014        beta=matpar->GetBeta();
     3015        meltingpoint=matpar->GetMeltingPoint();
     3016        this->parameters->FindParam(&dt,DtEnum);
     3017        Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input);
     3018
     3019        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
     3020        gauss=new GaussPenta(0,1,2,2);
     3021        for(ig=gauss->begin();ig<gauss->end();ig++){
     3022
     3023                gauss->GaussPoint(ig);
     3024
     3025                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
     3026                GetNodalFunctionsP1(&l1l6[0], gauss);
     3027
     3028                pressure_input->GetParameterValue(&pressure,gauss);
     3029                t_pmp=meltingpoint-beta*pressure;
     3030
     3031                scalar_ocean=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(t_pmp)/(heatcapacity*rho_ice);
     3032                if(dt) scalar_ocean=dt*scalar_ocean;
     3033
     3034                for(i=0;i<numdof;i++) pe->values[i]+=scalar_ocean*l1l6[i];
     3035        }
    29933036
    29943037        /*Clean up and return*/
     3038        delete gauss;
    29953039        return pe;
    29963040}
    … …  
    29993043ElementVector* Penta::CreatePVectorThermalSheet(void){
    30003044
     3045        /*Constants*/
     3046        const int  numdof=NUMVERTICES*NDOF1;
     3047
     3048        /*Intermediaries */
     3049        int        i,j,ig;
     3050        int        analysis_type,drag_type;
     3051        double     xyz_list[NUMVERTICES][3];
     3052        double     xyz_list_tria[NUMVERTICES2D][3];
     3053        double     Jdet2d,dt;
     3054        double     rho_ice,heatcapacity,geothermalflux_value;
     3055        double     basalfriction,alpha2,vx,vy,pressure;
     3056        double     pressure_list[3];
     3057        double     scalar;
     3058        double     l1l6[NUMVERTICES];
     3059        Friction*  friction=NULL;
     3060        GaussPenta* gauss=NULL;
     3061
    30013062        /* Geothermal flux on ice sheet base and basal friction */
    30023063        if (!IsOnBed() || IsOnShelf()) return NULL;
    30033064
    3004         /*Call Tria function*/
    3005         Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
    3006         ElementVector* pe=tria->CreatePVectorThermalSheet();
    3007         delete tria->matice; delete tria;
     3065        /*Initialize Element vector*/
     3066        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     3067
     3068        /*Retrieve all inputs and parameters*/
     3069        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     3070        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
     3071        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     3072        rho_ice=matpar->GetRhoIce();
     3073        heatcapacity=matpar->GetHeatCapacity();
     3074        this->parameters->FindParam(&dt,DtEnum);
     3075        Input* vx_input=inputs->GetInput(VxEnum);                         _assert_(vx_input);
     3076        Input* vy_input=inputs->GetInput(VyEnum);                         _assert_(vy_input);
     3077        Input* geothermalflux_input=inputs->GetInput(GeothermalFluxEnum); _assert_(geothermalflux_input);
     3078
     3079        /*Build frictoin element, needed later: */
     3080        inputs->GetParameterValue(&drag_type,DragTypeEnum);
     3081        if (drag_type!=2)_error_(" non-viscous friction not supported yet!");
     3082        friction=new Friction("3d",inputs,matpar,analysis_type);
     3083
     3084        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
     3085        gauss=new GaussPenta(0,1,2,2);
     3086        for(ig=gauss->begin();ig<gauss->end();ig++){
     3087
     3088                gauss->GaussPoint(ig);
     3089
     3090                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
     3091                GetNodalFunctionsP1(&l1l6[0], gauss);
     3092
     3093                geothermalflux_input->GetParameterValue(&geothermalflux_value,gauss);
     3094                friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
     3095                vx_input->GetParameterValue(&vx,gauss);
     3096                vy_input->GetParameterValue(&vy,gauss);
     3097                basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0));
     3098               
     3099                scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice);
     3100                if(dt) scalar=dt*scalar;
     3101
     3102                for(i=0;i<numdof;i++) pe->values[i]+=scalar*l1l6[i];
     3103        }
    30083104
    30093105        /*Clean up and return*/
     3106        delete gauss;
     3107        delete friction;
    30103108        return pe;
    30113109}

  • issm/trunk/src/c/objects/Elements/Tria.cpp

    r8409 r8410  
    25142514        /*Clean up and return*/
    25152515        delete gauss;
    2516         return pe;
    2517 }
    2518 /*}}}*/
    2519 /*FUNCTION Tria::CreatePVectorThermalShelf {{{1*/
    2520 ElementVector* Tria::CreatePVectorThermalShelf(void){
    2521        
    2522         /*Constants*/
    2523         const int  numdof=NUMVERTICES*NDOF1;
    2524 
    2525         /*Intermediaries */
    2526         int        i,ig;
    2527         double     Jdet;
    2528         double     mixed_layer_capacity,thermal_exchange_velocity;
    2529         double     rho_ice,rho_water,pressure,dt,scalar_ocean;
    2530         double     meltingpoint,beta,heatcapacity,t_pmp;
    2531         double     xyz_list[NUMVERTICES][3];
    2532         double     l1l2l3[NUMVERTICES];
    2533         GaussTria* gauss=NULL;
    2534 
    2535         /*Initialize Element vector*/
    2536         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
    2537 
    2538         /*Retrieve all inputs and parameters*/
    2539         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    2540         mixed_layer_capacity=matpar->GetMixedLayerCapacity();
    2541         thermal_exchange_velocity=matpar->GetThermalExchangeVelocity();
    2542         rho_water=matpar->GetRhoWater();
    2543         rho_ice=matpar->GetRhoIce();
    2544         heatcapacity=matpar->GetHeatCapacity();
    2545         beta=matpar->GetBeta();
    2546         meltingpoint=matpar->GetMeltingPoint();
    2547         this->parameters->FindParam(&dt,DtEnum);
    2548         Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input);
    2549 
    2550         /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    2551         gauss=new GaussTria(2);
    2552         for(ig=gauss->begin();ig<gauss->end();ig++){
    2553 
    2554                 gauss->GaussPoint(ig);
    2555 
    2556                 GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0], gauss);
    2557                 GetNodalFunctions(&l1l2l3[0], gauss);
    2558 
    2559                 pressure_input->GetParameterValue(&pressure,gauss);
    2560                 t_pmp=meltingpoint-beta*pressure;
    2561 
    2562                 scalar_ocean=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(t_pmp)/(heatcapacity*rho_ice);
    2563                 if(dt) scalar_ocean=dt*scalar_ocean;
    2564 
    2565                 for(i=0;i<numdof;i++) pe->values[i]+=scalar_ocean*l1l2l3[i];
    2566         }
    2567 
    2568         /*Clean up and return*/
    2569         delete gauss;
    2570         return pe;
    2571 }
    2572 /*}}}*/
    2573 /*FUNCTION Tria::CreatePVectorThermalSheet {{{1*/
    2574 ElementVector* Tria::CreatePVectorThermalSheet(void){
    2575 
    2576         /*Constants*/
    2577         const int  numdof=NUMVERTICES*NDOF1;
    2578 
    2579         /*Intermediaries */
    2580         int        i,ig;
    2581         int        analysis_type,drag_type;
    2582         double     xyz_list[NUMVERTICES][3];
    2583         double     Jdet,dt;
    2584         double     rho_ice,heatcapacity,geothermalflux_value;
    2585         double     basalfriction,alpha2,vx,vy,pressure;
    2586         double     pressure_list[3];
    2587         double     scalar;
    2588         double     l1l2l3[NUMVERTICES];
    2589         Friction*  friction=NULL;
    2590         GaussTria* gauss=NULL;
    2591 
    2592         /*Initialize Element vector*/
    2593         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
    2594 
    2595         /*Retrieve all inputs and parameters*/
    2596         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    2597         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    2598         rho_ice=matpar->GetRhoIce();
    2599         heatcapacity=matpar->GetHeatCapacity();
    2600         this->parameters->FindParam(&dt,DtEnum);
    2601         Input* vx_input=inputs->GetInput(VxEnum);                         _assert_(vx_input);
    2602         Input* vy_input=inputs->GetInput(VyEnum);                         _assert_(vy_input);
    2603         Input* vz_input=inputs->GetInput(VzEnum);                         _assert_(vz_input);
    2604         Input* geothermalflux_input=inputs->GetInput(GeothermalFluxEnum); _assert_(geothermalflux_input);
    2605 
    2606         /*Build frictoin element, needed later: */
    2607         inputs->GetParameterValue(&drag_type,DragTypeEnum);
    2608         if (drag_type!=2)_error_(" non-viscous friction not supported yet!");
    2609         friction=new Friction("3d",inputs,matpar,analysis_type);
    2610 
    2611         /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    2612         gauss=new GaussTria(2);
    2613         for(ig=gauss->begin();ig<gauss->end();ig++){
    2614 
    2615                 gauss->GaussPoint(ig);
    2616 
    2617                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0], gauss);
    2618                 GetNodalFunctions(&l1l2l3[0], gauss);
    2619 
    2620                 geothermalflux_input->GetParameterValue(&geothermalflux_value,gauss);
    2621                 friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
    2622                 vx_input->GetParameterValue(&vx,gauss);
    2623                 vy_input->GetParameterValue(&vy,gauss);
    2624                 basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0));
    2625                
    2626                 scalar=gauss->weight*Jdet*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice);
    2627                 if(dt) scalar=dt*scalar;
    2628 
    2629                 for(i=0;i<numdof;i++) pe->values[i]+=scalar*l1l2l3[i];
    2630         }
    2631 
    2632         /*Clean up and return*/
    2633         delete gauss;
    2634         delete friction;
    26352516        return pe;
    26362517}

  • issm/trunk/src/c/objects/Elements/Tria.h

    r8409 r8410  
    175175                ElementVector* CreatePVectorPrognostic_DG(void);
    176176                ElementVector* CreatePVectorSlope(void);
    177                 ElementVector* CreatePVectorThermalSheet(void);
    178                 ElementVector* CreatePVectorThermalShelf(void);
    179177                double  GetArea(void);
    180178                int     GetElementType(void);
Note: See TracChangeset for help on using the changeset viewer.