Changeset 8410 for issm/trunk/src/c/objects/Elements/Tria.cpp

Timestamp:

05/24/11 11:29:04 (14 years ago)

Author:

seroussi

Message:

moved thermal PVectors from Tria to Penta

File:

• issm/trunk/src/c/objects/Elements/Tria.cpp (modified) (1 diff)

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

@@ -2514 +2514 @@
         /*Clean up and return*/
         delete gauss;
-        return pe;
-}
-/*}}}*/
-/*FUNCTION Tria::CreatePVectorThermalShelf {{{1*/
-ElementVector* Tria::CreatePVectorThermalShelf(void){
-        
-        /*Constants*/
-        const int  numdof=NUMVERTICES*NDOF1;
-
-        /*Intermediaries */
-        int        i,ig;
-        double     Jdet;
-        double     mixed_layer_capacity,thermal_exchange_velocity;
-        double     rho_ice,rho_water,pressure,dt,scalar_ocean;
-        double     meltingpoint,beta,heatcapacity,t_pmp;
-        double     xyz_list[NUMVERTICES][3];
-        double     l1l2l3[NUMVERTICES];
-        GaussTria* gauss=NULL;
-
-        /*Initialize Element vector*/
-        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
-
-        /*Retrieve all inputs and parameters*/
-        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
-        mixed_layer_capacity=matpar->GetMixedLayerCapacity();
-        thermal_exchange_velocity=matpar->GetThermalExchangeVelocity();
-        rho_water=matpar->GetRhoWater();
-        rho_ice=matpar->GetRhoIce();
-        heatcapacity=matpar->GetHeatCapacity();
-        beta=matpar->GetBeta();
-        meltingpoint=matpar->GetMeltingPoint();
-        this->parameters->FindParam(&dt,DtEnum);
-        Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input);
-
-        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
-        gauss=new GaussTria(2);
-        for(ig=gauss->begin();ig<gauss->end();ig++){
-
-                gauss->GaussPoint(ig);
-
-                GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0], gauss);
-                GetNodalFunctions(&l1l2l3[0], gauss);
-
-                pressure_input->GetParameterValue(&pressure,gauss);
-                t_pmp=meltingpoint-beta*pressure;
-
-                scalar_ocean=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(t_pmp)/(heatcapacity*rho_ice);
-                if(dt) scalar_ocean=dt*scalar_ocean;
-
-                for(i=0;i<numdof;i++) pe->values[i]+=scalar_ocean*l1l2l3[i];
-        }
-
-        /*Clean up and return*/
-        delete gauss;
-        return pe;
-}
-/*}}}*/
-/*FUNCTION Tria::CreatePVectorThermalSheet {{{1*/
-ElementVector* Tria::CreatePVectorThermalSheet(void){
-
-        /*Constants*/
-        const int  numdof=NUMVERTICES*NDOF1;
-
-        /*Intermediaries */
-        int        i,ig;
-        int        analysis_type,drag_type;
-        double     xyz_list[NUMVERTICES][3];
-        double     Jdet,dt;
-        double     rho_ice,heatcapacity,geothermalflux_value;
-        double     basalfriction,alpha2,vx,vy,pressure;
-        double     pressure_list[3];
-        double     scalar;
-        double     l1l2l3[NUMVERTICES];
-        Friction*  friction=NULL;
-        GaussTria* gauss=NULL;
-
-        /*Initialize Element vector*/
-        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
-
-        /*Retrieve all inputs and parameters*/
-        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
-        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-        rho_ice=matpar->GetRhoIce();
-        heatcapacity=matpar->GetHeatCapacity();
-        this->parameters->FindParam(&dt,DtEnum);
-        Input* vx_input=inputs->GetInput(VxEnum);                         _assert_(vx_input);
-        Input* vy_input=inputs->GetInput(VyEnum);                         _assert_(vy_input);
-        Input* vz_input=inputs->GetInput(VzEnum);                         _assert_(vz_input);
-        Input* geothermalflux_input=inputs->GetInput(GeothermalFluxEnum); _assert_(geothermalflux_input);
-
-        /*Build frictoin element, needed later: */
-        inputs->GetParameterValue(&drag_type,DragTypeEnum);
-        if (drag_type!=2)_error_(" non-viscous friction not supported yet!");
-        friction=new Friction("3d",inputs,matpar,analysis_type);
-
-        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
-        gauss=new GaussTria(2);
-        for(ig=gauss->begin();ig<gauss->end();ig++){
-
-                gauss->GaussPoint(ig);
-
-                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0], gauss);
-                GetNodalFunctions(&l1l2l3[0], gauss);
-
-                geothermalflux_input->GetParameterValue(&geothermalflux_value,gauss);
-                friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
-                vx_input->GetParameterValue(&vx,gauss);
-                vy_input->GetParameterValue(&vy,gauss);
-                basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0));
-                
-                scalar=gauss->weight*Jdet*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice);
-                if(dt) scalar=dt*scalar;
-
-                for(i=0;i<numdof;i++) pe->values[i]+=scalar*l1l2l3[i];
-        }
-
-        /*Clean up and return*/
-        delete gauss;
-        delete friction;
         return pe;
 }