Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 8409) +++ /issm/trunk/src/c/Makefile.am (revision 8410) @@ -635,5 +635,5 @@ ./modules/ConstraintsStatex/RiftConstraintsState.cpp\ ./modules/ConstraintsStatex/ThermalConstraintsState.cpp\ - ./modules/ConstraintsStatex/MeltingIsPresent.cpp\ + ./modules/ConstraintsStatex/ThermalIsPresent.cpp\ ./modules/Responsex/Responsex.h\ ./modules/Responsex/Responsex.cpp\ @@ -1270,5 +1270,5 @@ ./modules/ConstraintsStatex/RiftConstraintsState.cpp\ ./modules/ConstraintsStatex/ThermalConstraintsState.cpp\ - ./modules/ConstraintsStatex/MeltingIsPresent.cpp\ + ./modules/ConstraintsStatex/ThermalIsPresent.cpp\ ./modules/Responsex/Responsex.h\ ./modules/Responsex/Responsex.cpp\ Index: /issm/trunk/src/c/modules/ConstraintsStatex/ConstraintsStateLocal.h =================================================================== --- /issm/trunk/src/c/modules/ConstraintsStatex/ConstraintsStateLocal.h (revision 8409) +++ /issm/trunk/src/c/modules/ConstraintsStatex/ConstraintsStateLocal.h (revision 8410) @@ -11,5 +11,5 @@ /*melting: */ void ThermalConstraintsState(Loads* loads, int* pconverged, int* pnum_unstable_constraints,int analysis_type); -int MeltingIsPresent(Loads* loads,int analysis_type); +int ThermalIsPresent(Loads* loads,int analysis_type); /*rifts module: */ Index: /issm/trunk/src/c/modules/ConstraintsStatex/ConstraintsStatex.cpp =================================================================== --- /issm/trunk/src/c/modules/ConstraintsStatex/ConstraintsStatex.cpp (revision 8409) +++ /issm/trunk/src/c/modules/ConstraintsStatex/ConstraintsStatex.cpp (revision 8410) @@ -35,5 +35,5 @@ RiftConstraintsState(&converged,&num_unstable_constraints,loads,min_mechanical_constraints,analysis_type); } - else if(MeltingIsPresent(loads,analysis_type)){ + else if(ThermalIsPresent(loads,analysis_type)){ ThermalConstraintsState(loads,&converged,&num_unstable_constraints,analysis_type); } Index: sm/trunk/src/c/modules/ConstraintsStatex/MeltingIsPresent.cpp =================================================================== --- /issm/trunk/src/c/modules/ConstraintsStatex/MeltingIsPresent.cpp (revision 8409) +++ (revision ) @@ -1,37 +1,0 @@ -/*!\file: MeltingIsPresent.cpp - * \brief - */ - -#ifdef HAVE_CONFIG_H - #include "config.h" -#else -#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" -#endif - -#include "./ConstraintsStateLocal.h" - -int MeltingIsPresent(Loads* loads,int configuration_type){ - - int i; - int found=0; - int mpi_found=0; - - for(i=0;iSize();i++){ - Object* object=(Object*)loads->GetObjectByOffset(i); - Load* load=(Load*)object; - if(load->InAnalysis(configuration_type)){ - if (object->Enum()==PengridEnum){ - found=1; - break; - } - } - } - - #ifdef _PARALLEL_ - MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); - MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD); - found=mpi_found; - #endif - - return found; -} Index: /issm/trunk/src/c/modules/ConstraintsStatex/ThermalIsPresent.cpp =================================================================== --- /issm/trunk/src/c/modules/ConstraintsStatex/ThermalIsPresent.cpp (revision 8410) +++ /issm/trunk/src/c/modules/ConstraintsStatex/ThermalIsPresent.cpp (revision 8410) @@ -0,0 +1,37 @@ +/*!\file: ThermalIsPresent.cpp + * \brief + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "./ConstraintsStateLocal.h" + +int ThermalIsPresent(Loads* loads,int configuration_type){ + + int i; + int found=0; + int mpi_found=0; + + for(i=0;iSize();i++){ + Object* object=(Object*)loads->GetObjectByOffset(i); + Load* load=(Load*)object; + if(load->InAnalysis(configuration_type)){ + if (object->Enum()==PengridEnum){ + found=1; + break; + } + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD); + found=mpi_found; + #endif + + return found; +} Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 8409) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 8410) @@ -2984,13 +2984,57 @@ 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;iGetMixedLayerCapacity(); + 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();igend();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;ivalues[i]+=scalar_ocean*l1l6[i]; + } /*Clean up and return*/ + delete gauss; return pe; } @@ -2999,13 +3043,67 @@ 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;iFindParam(&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();igend();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;ivalues[i]+=scalar*l1l6[i]; + } /*Clean up and return*/ + delete gauss; + delete friction; return pe; } Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 8409) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 8410) @@ -2514,123 +2514,4 @@ /*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();igend();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;ivalues[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();igend();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;ivalues[i]+=scalar*l1l2l3[i]; - } - - /*Clean up and return*/ - delete gauss; - delete friction; return pe; } Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 8409) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 8410) @@ -175,6 +175,4 @@ ElementVector* CreatePVectorPrognostic_DG(void); ElementVector* CreatePVectorSlope(void); - ElementVector* CreatePVectorThermalSheet(void); - ElementVector* CreatePVectorThermalShelf(void); double GetArea(void); int GetElementType(void);