Index: /issm/trunk-jpl/src/c/analyses/MeltingAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/MeltingAnalysis.cpp (revision 16809) +++ /issm/trunk-jpl/src/c/analyses/MeltingAnalysis.cpp (revision 16810) @@ -76,5 +76,5 @@ }/*}}}*/ ElementVector* MeltingAnalysis::CreatePVector(Element* element){/*{{{*/ -_error_("not implemented yet"); + return NULL; }/*}}}*/ void MeltingAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 16809) +++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 16810) @@ -116,5 +116,147 @@ }/*}}}*/ ElementVector* ThermalAnalysis::CreatePVector(Element* element){/*{{{*/ -_error_("not implemented yet"); + + /*compute all load vectors for this element*/ + ElementVector* pe1=CreatePVectorVolume(element); + ElementVector* pe2=CreatePVectorSheet(element); + ElementVector* pe3=CreatePVectorShelf(element); + ElementVector* pe =new ElementVector(pe1,pe2,pe3); + + /*clean-up and return*/ + delete pe1; + delete pe2; + delete pe3; + return pe; +}/*}}}*/ +ElementVector* ThermalAnalysis::CreatePVectorVolume(Element* element){/*{{{*/ + + /*Intermediaries*/ + int stabilization; + IssmDouble Jdet,phi,dt; + IssmDouble temperature; + IssmDouble tau_parameter,diameter; + IssmDouble u,v,w; + IssmDouble scalar_def,scalar_transient; + IssmDouble* xyz_list = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element vector*/ + ElementVector* pe = element->NewElementVector(); + IssmDouble* basis = xNew(numnodes); + IssmDouble* dbasis = xNew(3*numnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum); + IssmDouble thermalconductivity = element->GetMaterialParameter(MaterialsThermalconductivityEnum); + IssmDouble kappa = thermalconductivity/(rho_ice*heatcapacity); + element->FindParam(&dt,TimesteppingTimeStepEnum); + element->FindParam(&stabilization,ThermalStabilizationEnum); + Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input=element->GetInput(VzEnum); _assert_(vz_input); + Input* temperature_input = NULL; + if(reCast(dt)){temperature_input = element->GetInput(TemperatureEnum); _assert_(temperature_input);} + if(stabilization==2) diameter=element->MinEdgeLength(xyz_list); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(3); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctions(basis,gauss); + element->ViscousHeating(&phi,xyz_list,gauss,vx_input,vy_input,vz_input); + + scalar_def=phi/(rho_ice*heatcapacity)*Jdet*gauss->weight; + if(reCast(dt)) scalar_def=scalar_def*dt; + + for(int i=0;ivalues[i]+=scalar_def*basis[i]; + + /* Build transient now */ + if(reCast(dt)){ + temperature_input->GetInputValue(&temperature, gauss); + scalar_transient=temperature*Jdet*gauss->weight; + for(int i=0;ivalues[i]+=scalar_transient*basis[i]; + } + + if(stabilization==2){ + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + + vx_input->GetInputValue(&u,gauss); + vy_input->GetInputValue(&v,gauss); + vz_input->GetInputValue(&w,gauss); + + tau_parameter=element->StabilizationParameter(u,v,w,diameter,kappa); + + for(int i=0;ivalues[i]+=tau_parameter*scalar_def*(u*dbasis[0*3+i]+v*dbasis[1*3+i]+w*dbasis[2*3+i]); + if(reCast(dt)){ + for(int i=0;ivalues[i]+=tau_parameter*scalar_transient*(u*dbasis[0*3+i]+v*dbasis[1*3+i]+w*dbasis[2*3+i]); + } + } + } + + /*Clean up and return*/ + xDelete(basis); + xDelete(dbasis); + xDelete(xyz_list); + delete gauss; + return pe; + +}/*}}}*/ +ElementVector* ThermalAnalysis::CreatePVectorSheet(Element* element){/*{{{*/ + return NULL; +}/*}}}*/ +ElementVector* ThermalAnalysis::CreatePVectorShelf(Element* element){/*{{{*/ + + IssmDouble t_pmp,dt,Jdet,scalar_ocean,pressure; + IssmDouble *xyz_list_base = NULL; + + /*Get basal element*/ + if(!element->IsOnBed() || !element->IsFloating()) return NULL; + + /*Fetch number of nodes for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize vectors*/ + ElementVector* pe = element->NewElementVector(); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinatesBase(&xyz_list_base); + element->FindParam(&dt,TimesteppingTimeStepEnum); + Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input); + IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); + IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum); + IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum); + IssmDouble mixed_layer_capacity= element->GetMaterialParameter(MaterialsMixedLayerCapacityEnum); + IssmDouble thermal_exchange_vel= element->GetMaterialParameter(MaterialsThermalExchangeVelocityEnum); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGaussBase(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); + element->NodalFunctions(basis,gauss); + + pressure_input->GetInputValue(&pressure,gauss); + t_pmp=element->TMeltingPoint(pressure); + + scalar_ocean=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel*(t_pmp)/(heatcapacity*rho_ice); + if(reCast(dt)) scalar_ocean=dt*scalar_ocean; + + for(int i=0;ivalues[i]+=scalar_ocean*basis[i]; + } + + /*Clean up and return*/ + delete gauss; + xDelete(basis); + xDelete(xyz_list_base); + return pe; }/*}}}*/ void ThermalAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.h (revision 16809) +++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.h (revision 16810) @@ -23,4 +23,7 @@ ElementMatrix* CreateKMatrix(Element* element); ElementVector* CreatePVector(Element* element); + ElementVector* CreatePVectorVolume(Element* element); + ElementVector* CreatePVectorSheet(Element* element); + ElementVector* CreatePVectorShelf(Element* element); void GetSolutionFromInputs(Vector* solution,Element* element); void InputUpdateFromSolution(IssmDouble* solution,Element* element); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16809) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16810) @@ -53,5 +53,7 @@ virtual int FiniteElement(void)=0; virtual IssmDouble GetMaterialParameter(int enum_in)=0; + virtual IssmDouble MinEdgeLength(IssmDouble* xyz_list)=0; virtual void NodalFunctions(IssmDouble* basis,Gauss* gauss)=0; + virtual void NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss)=0; virtual void NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss)=0; virtual void NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss)=0; @@ -127,4 +129,6 @@ virtual void ReduceMatrices(ElementMatrix* Ke,ElementVector* pe)=0; virtual void ResetCoordinateSystem()=0; + virtual IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa)=0; + virtual void ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0; virtual int VelocityInterpolation()=0; virtual int PressureInterpolation()=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16809) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16810) @@ -1528,6 +1528,6 @@ } /*}}}*/ -/*FUNCTION Penta::GetStabilizationParameter {{{*/ -IssmDouble Penta::GetStabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){ +/*FUNCTION Penta::StabilizationParameter {{{*/ +IssmDouble Penta::StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){ /*Compute stabilization parameter*/ /*kappa=thermalconductivity/(rho_ice*hearcapacity) for thermal model*/ @@ -2465,5 +2465,5 @@ /*}}}*/ /*FUNCTION Penta::MinEdgeLength{{{*/ -IssmDouble Penta::MinEdgeLength(IssmDouble xyz_list[6][3]){ +IssmDouble Penta::MinEdgeLength(IssmDouble* xyz_list){ /*Return the minimum lenght of the nine egdes of the penta*/ @@ -2479,5 +2479,5 @@ /*Compute the length of this edge and compare it to the minimal length*/ - length=sqrt(pow(xyz_list[node0][0]-xyz_list[node1][0],2)+pow(xyz_list[node0][1]-xyz_list[node1][1],2)+pow(xyz_list[node0][2]-xyz_list[node1][2],2)); + length=sqrt(pow(xyz_list[node0*3+0]-xyz_list[node1*3+0],2)+pow(xyz_list[node0*3+1]-xyz_list[node1*3+1],2)+pow(xyz_list[node0*3+2]-xyz_list[node1*3+2],2)); if(lengthEnum()==GaussPentaEnum); this->GetNodalFunctions(basis,(GaussPenta*)gauss); + +} +/*}}}*/ +/*FUNCTION Penta::NodalFunctionsDerivatives{{{*/ +void Penta::NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){ + + _assert_(gauss->Enum()==GaussPentaEnum); + this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss); } @@ -3379,4 +3387,21 @@ } /*}}}*/ +/*FUNCTION Penta::ViscousHeating{{{*/ +void Penta::ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){ + + /*Intermediaries*/ + IssmDouble phi; + IssmDouble viscosity; + IssmDouble epsilon[6]; + + _assert_(gauss->Enum()==GaussPentaEnum); + this->GetStrainRate3d(&epsilon[0],xyz_list,(GaussPenta*)gauss,vx_input,vy_input,vz_input); + material->GetViscosity3dFS(&viscosity,&epsilon[0]); + GetPhi(&phi,&epsilon[0],viscosity); + + /*Assign output pointer*/ + *pphi = phi; +} +/*}}}*/ /*FUNCTION Penta::VelocityInterpolation{{{*/ int Penta::VelocityInterpolation(void){ @@ -3742,5 +3767,5 @@ Input* vym_input=inputs->GetInput(VyMeshEnum); _assert_(vym_input); Input* vzm_input=inputs->GetInput(VzMeshEnum); _assert_(vzm_input); - if (stabilization==2) diameter=MinEdgeLength(xyz_list); + if (stabilization==2) diameter=MinEdgeLength(&xyz_list[0][0]); /* Start looping on the number of gaussian points: */ @@ -3825,5 +3850,5 @@ else if(stabilization==2){ GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss); - tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,kappa/rho_ice); + tau_parameter=StabilizationParameter(u-um,v-vm,w-wm,diameter,kappa/rho_ice); for(i=0;iGetInput(VyMeshEnum); _assert_(vym_input); Input* vzm_input=inputs->GetInput(VzMeshEnum); _assert_(vzm_input); - if (stabilization==2) diameter=MinEdgeLength(xyz_list); + if (stabilization==2) diameter=MinEdgeLength(&xyz_list[0][0]); /* Start looping on the number of gaussian points: */ @@ -4054,5 +4079,5 @@ else if(stabilization==2){ GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss); - tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,kappa); + tau_parameter=StabilizationParameter(u-um,v-vm,w-wm,diameter,kappa); for(i=0;iGetEnthalpyDiffusionParameterVolume(enthalpypicard,pressure_p); kappa/=rho_ice; - tau_parameter=GetStabilizationParameter(u,v,w,diameter,kappa); + tau_parameter=StabilizationParameter(u,v,w,diameter,kappa); for(i=0;ivalues[i]+=tau_parameter*scalar_def*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]); @@ -4459,5 +4484,5 @@ Input* temperature_input=NULL; if(reCast(dt)) temperature_input=inputs->GetInput(TemperatureEnum); _assert_(inputs); - if(stabilization==2) diameter=MinEdgeLength(xyz_list); + if(stabilization==2) diameter=MinEdgeLength(&xyz_list[0][0]); /* Start looping on the number of gaussian points: */ @@ -4493,5 +4518,5 @@ vz_input->GetInputValue(&w, gauss); - tau_parameter=GetStabilizationParameter(u,v,w,diameter,kappa); + tau_parameter=StabilizationParameter(u,v,w,diameter,kappa); for(i=0;ivalues[i]+=tau_parameter*scalar_def*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]); @@ -7615,5 +7640,5 @@ /*Find minimal length and B*/ rigidity=material->GetB(); - diameter=MinEdgeLength(xyz_list); + diameter=MinEdgeLength(&xyz_list[0][0]); gauss=new GaussPenta(); @@ -8756,5 +8781,5 @@ /*Find minimal length*/ - diameter=MinEdgeLength(xyz_list); + diameter=MinEdgeLength(&xyz_list[0][0]); gauss=new GaussPenta(); Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16809) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16810) @@ -124,4 +124,5 @@ IssmDouble TimeAdapt(); void ViscousHeatingCreateInput(void); + void ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); #ifdef _HAVE_RESPONSES_ @@ -224,5 +225,4 @@ void GetPhi(IssmDouble* phi, IssmDouble* epsilon, IssmDouble viscosity); void GetQuadNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]); - IssmDouble GetStabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa); void GetStrainRate3dHO(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input); void GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input); @@ -247,10 +247,12 @@ ElementVector* NewElementVector(int approximation_enum); void NodalFunctions(IssmDouble* basis,Gauss* gauss); + void NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss); void NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss); void NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss); - IssmDouble MinEdgeLength(IssmDouble xyz_list[6][3]); + IssmDouble MinEdgeLength(IssmDouble* xyz_list); void SetClone(int* minranks); Tria* SpawnTria(int location); void SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]); + IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa); void TransformLoadVectorCoord(ElementVector* pe,int transformenum); void TransformLoadVectorCoord(ElementVector* pe,int* transformenum_list); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16809) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16810) @@ -110,7 +110,9 @@ void JacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; void JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");}; + IssmDouble MinEdgeLength(IssmDouble* xyz_list){_error_("not implemented yet");}; void NodalFunctions(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");}; void NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");}; void NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");}; + void NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; bool NoIceInElement(){_error_("not implemented yet");}; void NormalBase(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");}; @@ -118,4 +120,5 @@ int NumberofNodesPressure(void){_error_("not implemented yet");}; Element* SpawnBasalElement(void){_error_("not implemented yet");}; + IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");}; int PressureInterpolation(void){_error_("not implemented yet");}; int VelocityInterpolation(void){_error_("not implemented yet");}; @@ -132,4 +135,5 @@ Gauss* NewGaussBase(int order){_error_("not implemented yet");}; ElementVector* NewElementVector(int approximation_enum){_error_("not implemented yet");}; + void ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented yet");}; #ifdef _HAVE_THERMAL_ void UpdateBasalConstraintsEnthalpy(void){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16809) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16810) @@ -2213,4 +2213,12 @@ _assert_(gauss->Enum()==GaussTriaEnum); this->GetNodalFunctions(basis,(GaussTria*)gauss); + +} +/*}}}*/ +/*FUNCTION Tria::NodalFunctionsDerivatives{{{*/ +void Tria::NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){ + + _assert_(gauss->Enum()==GaussTriaEnum); + this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss); } Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16809) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16810) @@ -269,4 +269,5 @@ void JacobianDeterminant(IssmDouble* pJdet, IssmDouble* xyz_list,Gauss* gauss); void JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");}; + IssmDouble MinEdgeLength(IssmDouble* xyz_list){_error_("not implemented yet");}; Gauss* NewGauss(void); Gauss* NewGauss(int order); @@ -274,8 +275,10 @@ ElementVector* NewElementVector(int approximation_enum); void NodalFunctions(IssmDouble* basis,Gauss* gauss); + void NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss); void NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss); void NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss); void SetClone(int* minranks); Seg* SpawnSeg(int index1,int index2); + IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");}; void SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]); void TransformLoadVectorCoord(ElementVector* pe,int transformenum); @@ -287,4 +290,5 @@ void TransformSolutionCoord(IssmDouble* values,int numnodes,int transformenum){_error_("not implemented yet");}; void TransformSolutionCoord(IssmDouble* values,int numnodes,int* transformenum_list){_error_("not implemented yet");}; + void ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented yet");}; #ifdef _HAVE_STRESSBALANCE_ ElementMatrix* CreateKMatrixStressbalanceSSA(void);