Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5931) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5932) @@ -789,8 +789,12 @@ break; case ThermalAnalysisEnum: - CreatePVectorThermal( pg); + pe=CreatePVectorThermal(); + if(pe) pe->AddToGlobal(pg,pf); + if(pe) delete pe; break; case MeltingAnalysisEnum: - CreatePVectorMelting( pg); + pe=CreatePVectorMelting(); + if(pe) pe->AddToGlobal(pg,pf); + if(pe) delete pe; break; default: @@ -3736,6 +3740,6 @@ /*}}}*/ /*FUNCTION Penta::CreatePVectorMelting {{{1*/ -void Penta::CreatePVectorMelting( Vec pg){ - return; +ElementVector* Penta::CreatePVectorMelting(void){ + return NULL; } /*}}}*/ @@ -3777,32 +3781,34 @@ /*}}}*/ /*FUNCTION Penta::CreatePVectorThermal {{{1*/ -void Penta::CreatePVectorThermal( Vec pg){ - - /*indexing: */ +ElementVector* Penta::CreatePVectorThermal(void){ + + /*compute all load vectorsfor this element*/ + ElementVector* pe1=CreatePVectorThermalVolume(); + ElementVector* pe2=CreatePVectorThermalSheet(); + ElementVector* pe3=CreatePVectorThermalShelf(); + ElementVector* pe =new ElementVector(pe1,pe2,pe3); + + /*clean-up and return*/ + delete pe1; + delete pe2; + delete pe3; + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorThermalVolume {{{1*/ +ElementVector* Penta::CreatePVectorThermalVolume(void){ + + const int numdof=NUMVERTICES*NDOF1; int i,j; + int ig; int found=0; - - const int numdof=NUMVERTICES*NDOF1; - int* doflist=NULL; - - /*Grid data: */ double xyz_list[NUMVERTICES][3]; - - /* gaussian points: */ - int ig; GaussPenta *gauss=NULL; - double temperature_list[NUMVERTICES]; double temperature; - - /*Material properties: */ double gravity,rho_ice,rho_water; double mixed_layer_capacity,heatcapacity; double beta,meltingpoint,thermal_exchange_velocity; - - /* element parameters: */ int friction_type; - - /*matrices: */ double P_terms[numdof]={0.0}; double L[numdof]; @@ -3812,5 +3818,4 @@ double epsilon_sqr[3][3]; double epsilon_matrix[3][3]; - double Jdet; double viscosity; @@ -3822,20 +3827,14 @@ double scalar_ocean; double scalar_transient; - - /*Collapsed formulation: */ Tria* tria=NULL; double dt; - /*retrieve some parameters: */ + /*Initialize Element vector and return if necessary*/ + if(IsOnWater()) return NULL; + ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); this->parameters->FindParam(&dt,DtEnum); - - /*If on water, skip: */ - if(IsOnWater())return; - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - - /*recovre material parameters: */ rho_water=matpar->GetRhoWater(); rho_ice=matpar->GetRhoIce(); @@ -3844,6 +3843,4 @@ beta=matpar->GetBeta(); meltingpoint=matpar->GetMeltingPoint(); - - /*Retrieve all inputs we will be needing: */ Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); @@ -3860,22 +3857,14 @@ gauss->GaussPoint(ig); - /*Compute strain rate and viscosity: */ this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); - - /* Get Jacobian determinant: */ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - - /* Get nodal functions */ GetNodalFunctionsP1(&L[0], gauss); - - /*Build deformational heating: */ GetPhi(&phi, &epsilon[0], viscosity); - /*Build pe_gaussian */ scalar_def=phi/(rho_ice*heatcapacity)*Jdet*gauss->weight; if(dt) scalar_def=scalar_def*dt; - for(i=0;ivalues[i]+=scalar_def*L[i]; /* Build transient now */ @@ -3883,30 +3872,41 @@ temperature_input->GetParameterValue(&temperature, gauss); scalar_transient=temperature*Jdet*gauss->weight; - for(i=0;ivalues[i]+=scalar_transient*L[i]; + } + } + + /*Clean up and return*/ + delete gauss; + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorThermalShelf {{{1*/ +ElementVector* Penta::CreatePVectorThermalShelf(void){ /* Ice/ocean heat exchange flux on ice shelf base */ - if(IsOnBed() && IsOnShelf()){ - - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVectorThermalShelf(pg); - delete tria->matice; delete tria; - } + if (!IsOnBed() || !IsOnShelf() || IsOnWater()) return NULL; + + /*Call Tria function*/ + Tria* tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + ElementVector* pe=tria->CreatePVectorThermalShelf(); + delete tria->matice; delete tria; + + /*Clean up and return*/ + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorThermalSheet {{{1*/ +ElementVector* Penta::CreatePVectorThermalSheet(void){ /* Geothermal flux on ice sheet base and basal friction */ - if(IsOnBed() && !IsOnShelf()){ - - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVectorThermalSheet(pg); - delete tria->matice; delete tria; - } - - xfree((void**)&doflist); - delete gauss; - + if (!IsOnBed() || IsOnShelf() || IsOnWater()) return NULL; + + /*Call Tria function*/ + Tria* tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + ElementVector* pe=tria->CreatePVectorThermalSheet(); + delete tria->matice; delete tria; + + /*Clean up and return*/ + return pe; } /*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 5931) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 5932) @@ -169,8 +169,11 @@ ElementVector* CreatePVectorAdjointStokes(void); void CreatePVectorDiagnosticVert( Vec pg); - void CreatePVectorMelting( Vec pg); + ElementVector* CreatePVectorMelting(void); ElementVector* CreatePVectorPrognostic(void); ElementVector* CreatePVectorSlope(void); - void CreatePVectorThermal( Vec pg); + ElementVector* CreatePVectorThermal(void); + ElementVector* CreatePVectorThermalVolume(void); + ElementVector* CreatePVectorThermalShelf(void); + ElementVector* CreatePVectorThermalSheet(void); void GetDofList(int** pdoflist,int approximation_enum,int setenum); void GetDofList1(int* doflist); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5931) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5932) @@ -4388,12 +4388,9 @@ /*}}}*/ /*FUNCTION Tria::CreatePVectorThermalShelf {{{1*/ -void Tria::CreatePVectorThermalShelf( Vec pg){ - - int i; +ElementVector* Tria::CreatePVectorThermalShelf(void){ const int numdof=NUMVERTICES*NDOF1; - int* doflist=NULL; + int i; double xyz_list[NUMVERTICES][3]; - double mixed_layer_capacity; double thermal_exchange_velocity; @@ -4403,25 +4400,19 @@ double beta; double meltingpoint; - - /*inputs: */ double dt; double pressure; - - /* gaussian points: */ int ig; GaussTria* gauss=NULL; - - /*matrices: */ double Jdet; - double P_terms[numdof]={0.0}; double l1l2l3[NUMVERTICES]; double t_pmp; double scalar_ocean; - /* Get node coordinates and dof list: */ + /*Initialize Element vector and return if necessary*/ + if(IsOnWater()) return NULL; + ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters); + + /*Retrieve all inputs and parameters*/ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - - //recover material parameters mixed_layer_capacity=matpar->GetMixedLayerCapacity(); thermal_exchange_velocity=matpar->GetThermalExchangeVelocity(); @@ -4431,11 +4422,5 @@ beta=matpar->GetBeta(); meltingpoint=matpar->GetMeltingPoint(); - - /*retrieve some solution parameters: */ this->parameters->FindParam(&dt,DtEnum); - - /* Ice/ocean heat exchange flux on ice shelf base */ - - /*Retrieve all inputs we will be needing: */ Input* pressure_input=inputs->GetInput(PressureEnum); ISSMASSERT(pressure_input); @@ -4446,15 +4431,10 @@ gauss->GaussPoint(ig); - //Get the Jacobian determinant GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0], gauss); - - /*Get nodal functions values: */ GetNodalFunctions(&l1l2l3[0], gauss); - /*Get geothermal flux and basal friction */ pressure_input->GetParameterValue(&pressure,gauss); t_pmp=meltingpoint-beta*pressure; - /*Calculate scalar parameter*/ scalar_ocean=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(t_pmp)/(heatcapacity*rho_ice); if(dt){ @@ -4462,30 +4442,21 @@ } - for(i=0;i<3;i++){ - P_terms[i]+=scalar_ocean*l1l2l3[i]; - } - } - - /*Add pe_g to global vector pg: */ - VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES); + for(i=0;ivalues[i]+=scalar_ocean*l1l2l3[i]; + } /*Clean up and return*/ delete gauss; - xfree((void**)&doflist); + return pe; } /*}}}*/ /*FUNCTION Tria::CreatePVectorThermalSheet {{{1*/ -void Tria::CreatePVectorThermalSheet( Vec pg){ - +ElementVector* Tria::CreatePVectorThermalSheet(void){ + + const int numdof=NUMVERTICES*NDOF1; int i; - - const int numdof=NUMVERTICES*NDOF1; - int* doflist=NULL; + int ig; double xyz_list[NUMVERTICES][3]; - double rho_ice; double heatcapacity; - - /*inputs: */ double dt; double pressure_list[3]; @@ -4496,30 +4467,25 @@ double alpha2,vx,vy; double geothermalflux_value; - - /* gaussian points: */ - int ig; GaussTria* gauss=NULL; - - /*matrices: */ double Jdet; double P_terms[numdof]={0.0}; double l1l2l3[NUMVERTICES]; double scalar; - int analysis_type; - /*retrive parameters: */ + /*Initialize Element vector and return if necessary*/ + if(IsOnWater()) return NULL; + ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - - //recover material parameters rho_ice=matpar->GetRhoIce(); heatcapacity=matpar->GetHeatCapacity(); - - /*retrieve some parameters: */ this->parameters->FindParam(&dt,DtEnum); + Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); + Input* vz_input=inputs->GetInput(VzEnum); ISSMASSERT(vz_input); + Input* geothermalflux_input=inputs->GetInput(GeothermalFluxEnum); ISSMASSERT(geothermalflux_input); /*Build frictoin element, needed later: */ @@ -4528,12 +4494,4 @@ friction=new Friction("3d",inputs,matpar,analysis_type); - /*Retrieve all inputs we will be needing: */ - Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); - Input* vz_input=inputs->GetInput(VzEnum); ISSMASSERT(vz_input); - Input* geothermalflux_input=inputs->GetInput(GeothermalFluxEnum); ISSMASSERT(geothermalflux_input); - - /* Ice/ocean heat exchange flux on ice shelf base */ - /* Start looping on the number of gauss 2d (nodes on the bedrock) */ gauss=new GaussTria(2); @@ -4542,17 +4500,10 @@ gauss->GaussPoint(ig); - //Get the Jacobian determinant GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0], gauss); - - /*Get nodal functions values: */ GetNodalFunctions(&l1l2l3[0], gauss); - /*Get geothermal flux and basal friction */ geothermalflux_input->GetParameterValue(&geothermalflux_value,gauss); - - /*Friction: */ friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum); - /*Calculate scalar parameter*/ scalar=gauss->weight*Jdet*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice); if(dt){ @@ -4560,17 +4511,11 @@ } - for(i=0;i<3;i++){ - P_terms[i]+=scalar*l1l2l3[i]; - } - - } - - /*Add pe_g to global vector pg: */ - VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES); + for(i=0;ivalues[i]+=scalar*l1l2l3[i]; + } /*Clean up and return*/ delete gauss; delete friction; - xfree((void**)&doflist); + return pe; } /*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 5931) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 5932) @@ -152,6 +152,6 @@ ElementVector* CreatePVectorPrognostic_DG(void); ElementVector* CreatePVectorSlope(void); - void CreatePVectorThermalSheet( Vec pg); - void CreatePVectorThermalShelf( Vec pg); + ElementVector* CreatePVectorThermalSheet(void); + ElementVector* CreatePVectorThermalShelf(void); double GetArea(void); int GetElementType(void);