Changeset 8410 for issm/trunk/src/c/objects/Elements/Penta.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/Penta.cpp (modified) (2 diffs)

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

@@ -2984 +2984 @@
 ElementVector* Penta::CreatePVectorThermalShelf(void){
 
+        /*Constants*/
+        const int  numdof=NUMVERTICES*NDOF1;
+
+        /*Intermediaries */
+        int        i,j,ig;
+        double     Jdet2d;
+        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     xyz_list_tria[NUMVERTICES2D][3];
+        double     l1l6[NUMVERTICES];
+        GaussPenta* gauss=NULL;
+
+        /*Initialize Element vector*/
+        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
+
         /* Ice/ocean heat exchange flux on ice shelf base */
         if (!IsOnBed() || !IsOnShelf()) return NULL;
 
-        /*Call Tria function*/
-        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
-        ElementVector* pe=tria->CreatePVectorThermalShelf();
-        delete tria->matice; delete tria;
+        /*Retrieve all inputs and parameters*/
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
+        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 GaussPenta(0,1,2,2);
+        for(ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
+                GetNodalFunctionsP1(&l1l6[0], gauss);
+
+                pressure_input->GetParameterValue(&pressure,gauss);
+                t_pmp=meltingpoint-beta*pressure;
+
+                scalar_ocean=gauss->weight*Jdet2d*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*l1l6[i];
+        }
 
         /*Clean up and return*/
+        delete gauss;
         return pe;
 }
@@ -2999 +3043 @@
 ElementVector* Penta::CreatePVectorThermalSheet(void){
 
+        /*Constants*/
+        const int  numdof=NUMVERTICES*NDOF1;
+
+        /*Intermediaries */
+        int        i,j,ig;
+        int        analysis_type,drag_type;
+        double     xyz_list[NUMVERTICES][3];
+        double     xyz_list_tria[NUMVERTICES2D][3];
+        double     Jdet2d,dt;
+        double     rho_ice,heatcapacity,geothermalflux_value;
+        double     basalfriction,alpha2,vx,vy,pressure;
+        double     pressure_list[3];
+        double     scalar;
+        double     l1l6[NUMVERTICES];
+        Friction*  friction=NULL;
+        GaussPenta* gauss=NULL;
+
         /* Geothermal flux on ice sheet base and basal friction */
         if (!IsOnBed() || IsOnShelf()) return NULL;
 
-        /*Call Tria function*/
-        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
-        ElementVector* pe=tria->CreatePVectorThermalSheet();
-        delete tria->matice; delete tria;
+        /*Initialize Element vector*/
+        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
+
+        /*Retrieve all inputs and parameters*/
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
+        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* 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 GaussPenta(0,1,2,2);
+        for(ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
+                GetNodalFunctionsP1(&l1l6[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*Jdet2d*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice);
+                if(dt) scalar=dt*scalar;
+
+                for(i=0;i<numdof;i++) pe->values[i]+=scalar*l1l6[i];
+        }
 
         /*Clean up and return*/
+        delete gauss;
+        delete friction;
         return pe;
 }