Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16981) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16982) @@ -105,5 +105,4 @@ virtual void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum)=0; virtual ElementMatrix* CreateKMatrix(void)=0; - virtual ElementVector* CreatePVector(void)=0; virtual void CreateDVector(Vector* df)=0; virtual void CreateJacobianMatrix(Matrix* Jff)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16981) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16982) @@ -540,147 +540,4 @@ delete De; } -} -/*}}}*/ -/*FUNCTION Penta::CreatePVector(void) {{{*/ -ElementVector* Penta::CreatePVector(void){ - - /*retrieve parameters: */ - int analysis_type; - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*if debugging mode, check that all pointers exist {{{*/ - _assert_(this->nodes && this->material && this->matpar && this->verticalneighbors && this->parameters && this->inputs); - /*}}}*/ - - /*Skip if water element*/ - if(NoIceInElement()) return NULL; - - /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ - switch(analysis_type){ - #ifdef _HAVE_STRESSBALANCE_ - case StressbalanceAnalysisEnum: - return CreatePVectorStressbalanceHoriz(); - break; - case StressbalanceSIAAnalysisEnum: - return CreatePVectorStressbalanceSIA(); - break; - case StressbalanceVerticalAnalysisEnum: - return CreatePVectorStressbalanceVert(); - break; - #endif - #ifdef _HAVE_CONTROL_ - case AdjointHorizAnalysisEnum: - return CreatePVectorAdjointHoriz(); - break; - #endif - #ifdef _HAVE_THERMAL_ - case ThermalAnalysisEnum: - return CreatePVectorThermal(); - break; - case EnthalpyAnalysisEnum: - return CreatePVectorEnthalpy(); - break; - case MeltingAnalysisEnum: - return CreatePVectorMelting(); - break; - #endif - case L2ProjectionBaseAnalysisEnum: - return CreatePVectorL2ProjectionBase(); - break; - case MasstransportAnalysisEnum: - return CreatePVectorMasstransport(); - break; - case FreeSurfaceTopAnalysisEnum: - return CreatePVectorFreeSurfaceTop(); - break; - case FreeSurfaceBaseAnalysisEnum: - return CreatePVectorFreeSurfaceBase(); - break; - #ifdef _HAVE_BALANCED_ - case BalancethicknessAnalysisEnum: - return CreatePVectorBalancethickness(); - break; - #endif - #ifdef _HAVE_HYDROLOGY_ - case HydrologyDCInefficientAnalysisEnum: - return CreatePVectorHydrologyDCInefficient(); - break; - case HydrologyDCEfficientAnalysisEnum: - return CreatePVectorHydrologyDCEfficient(); - break; - case L2ProjectionEPLAnalysisEnum: - return CreatePVectorL2ProjectionEPL(); - break; - #endif - default: - _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet"); - } - - /*make compiler happy*/ - return NULL; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorMasstransport {{{*/ -ElementVector* Penta::CreatePVectorMasstransport(void){ - - if (!IsOnBed()) return NULL; - - /*Depth Averaging Vx and Vy*/ - this->InputDepthAverageAtBase(VxEnum,VxAverageEnum); - this->InputDepthAverageAtBase(VyEnum,VyAverageEnum); - - /*Call Tria function*/ - Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. - ElementVector* pe=tria->CreatePVectorMasstransport(); - delete tria->material; delete tria; - - /*Delete Vx and Vy averaged*/ - this->inputs->DeleteInput(VxAverageEnum); - this->inputs->DeleteInput(VyAverageEnum); - - /*Clean up and return*/ - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorFreeSurfaceTop {{{*/ -ElementVector* Penta::CreatePVectorFreeSurfaceTop(void){ - - if(!IsOnSurface()) return NULL; - - /*Call Tria function*/ - Tria* tria=(Tria*)SpawnTria(1); //lower face is 0, upper face is 1. - ElementVector* pe=tria->CreatePVectorFreeSurfaceTop(); - delete tria->material; delete tria; - - /*Clean up and return*/ - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorFreeSurfaceBase {{{*/ -ElementVector* Penta::CreatePVectorFreeSurfaceBase(void){ - - if(!IsOnBed()) return NULL; - - /*Call Tria function*/ - Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. - ElementVector* pe=tria->CreatePVectorFreeSurfaceBase(); - delete tria->material; delete tria; - - /*Clean up and return*/ - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorL2ProjectionBase {{{*/ -ElementVector* Penta::CreatePVectorL2ProjectionBase(void){ - - if (!IsOnBed()) return NULL; - - /*Call Tria function*/ - Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. - ElementVector* pe=tria->CreatePVectorL2Projection(); - delete tria->material; delete tria; - - /*clean up and return*/ - return pe; } /*}}}*/ @@ -4031,504 +3888,4 @@ } /*}}}*/ -/*FUNCTION Penta::CreatePVectorEnthalpy {{{*/ -ElementVector* Penta::CreatePVectorEnthalpy(void){ - - /*compute all load vectors for this element*/ - ElementVector* pe1=CreatePVectorEnthalpyVolume(); - ElementVector* pe2=CreatePVectorEnthalpySheet(); - ElementVector* pe3=CreatePVectorEnthalpyShelf(); - ElementVector* pe =new ElementVector(pe1,pe2,pe3); - - /*clean-up and return*/ - delete pe1; - delete pe2; - delete pe3; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorEnthalpyVolume {{{*/ -ElementVector* Penta::CreatePVectorEnthalpyVolume(void){ - - /*Constants*/ - const int numdof=NUMVERTICES*NDOF1; - - /*Intermediaries*/ - int i; - int stabilization; - IssmDouble Jdet,phi,dt; - IssmDouble rho_ice,heatcapacity; - IssmDouble thermalconductivity,kappa; - IssmDouble viscosity,pressure,enthalpy; - IssmDouble enthalpypicard[NUMVERTICES], pressure_p[NUMVERTICES]; - IssmDouble tau_parameter,diameter; - IssmDouble u,v,w; - IssmDouble scalar_def,scalar_transient; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble L[numdof]; - IssmDouble dbasis[3][6]; - IssmDouble epsilon[6]; - GaussPenta *gauss=NULL; - - /*Initialize Element vector*/ - ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - rho_ice=matpar->GetRhoIce(); - heatcapacity=matpar->GetHeatCapacity(); - thermalconductivity=matpar->GetThermalConductivity(); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - this->parameters->FindParam(&stabilization,ThermalStabilizationEnum); - 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* enthalpy_input=NULL; - if(reCast(dt)){ - enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input); - } - if (stabilization==2){ - diameter=MinEdgeLength(&xyz_list[0][0]); - } - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(3,3); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctionsP1(&L[0], gauss); - - this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - material->GetViscosity3dFS(&viscosity,&epsilon[0]); - GetPhi(&phi, &epsilon[0], viscosity); - - scalar_def=phi/rho_ice*Jdet*gauss->weight; - if(reCast(dt)) scalar_def=scalar_def*dt; - - /*TODO: add -beta*laplace T_m(p)*/ - for(i=0;ivalues[i]+=scalar_def*L[i]; - - /* Build transient now */ - if(reCast(dt)){ - enthalpy_input->GetInputValue(&enthalpy, gauss); - scalar_transient=enthalpy*Jdet*gauss->weight; - for(i=0;ivalues[i]+=scalar_transient*L[i]; - } - - if(stabilization==2){ - GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss); - - vx_input->GetInputValue(&u, gauss); - vy_input->GetInputValue(&v, gauss); - vz_input->GetInputValue(&w, gauss); - GetInputListOnVertices(&enthalpypicard[0],EnthalpyPicardEnum); - GetInputListOnVertices(&pressure_p[0],PressureEnum); - kappa=matpar->GetEnthalpyDiffusionParameterVolume(NUMVERTICES,enthalpypicard,pressure_p); - kappa/=rho_ice; - 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]); - if(reCast(dt)){ - for(i=0;ivalues[i]+=tau_parameter*scalar_transient*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]); - } - } - } - - /*Clean up and return*/ - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorEnthalpyShelf {{{*/ -ElementVector* Penta::CreatePVectorEnthalpyShelf(void){ - - /*Constants*/ - const int numdof=NUMVERTICES*NDOF1; - - /*Intermediaries */ - int i,j; - IssmDouble Jdet2d; - IssmDouble heatcapacity,h_pmp; - IssmDouble mixed_layer_capacity,thermal_exchange_velocity; - IssmDouble rho_ice,rho_water,pressure,dt,scalar_ocean; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_tria[NUMVERTICES2D][3]; - IssmDouble basis[NUMVERTICES]; - GaussPenta* gauss=NULL; - - /* Ice/ocean heat exchange flux on ice shelf base */ - if (!IsOnBed() || !IsFloating()) return NULL; - - /*Initialize Element vector*/ - ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - for(i=0;iGetMixedLayerCapacity(); - thermal_exchange_velocity=matpar->GetThermalExchangeVelocity(); - rho_water=matpar->GetRhoWater(); - rho_ice=matpar->GetRhoIce(); - heatcapacity=matpar->GetHeatCapacity(); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - 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(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); - GetNodalFunctionsP1(&basis[0], gauss); - - pressure_input->GetInputValue(&pressure,gauss); - h_pmp=matpar->PureIceEnthalpy(pressure); - - scalar_ocean=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(h_pmp)/(rho_ice*heatcapacity); - if(reCast(dt)) scalar_ocean=dt*scalar_ocean; - - for(i=0;ivalues[i]+=scalar_ocean*basis[i]; - } - - /*Clean up and return*/ - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorEnthalpySheet {{{*/ -ElementVector* Penta::CreatePVectorEnthalpySheet(void){ - - /*Constants*/ - const int numdof=NUMVERTICES*NDOF1; - - /*Intermediaries */ - int i,j; - int analysis_type; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_tria[NUMVERTICES2D][3]={0.0}; - IssmDouble Jdet2d,dt; - IssmDouble rho_ice,heatcapacity; - IssmDouble geothermalflux_value, heatflux_value; - IssmDouble basalfriction,alpha2,vx,vy,vz; - IssmDouble scalar,enthalpy,enthalpyup; - IssmDouble pressure,pressureup; - IssmDouble watercolumn; - IssmDouble basis[NUMVERTICES]; - Friction* friction=NULL; - GaussPenta* gauss=NULL; - GaussPenta* gaussup=NULL; - - /* Geothermal flux on ice sheet base and basal friction */ - if (!IsOnBed() || IsFloating()) return NULL; - - /*Initialize Element vector*/ - ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - for(i=0;iFindParam(&analysis_type,AnalysisTypeEnum); - rho_ice=matpar->GetRhoIce(); - heatcapacity=matpar->GetHeatCapacity(); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - 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* enthalpy_input=inputs->GetInput(EnthalpyPicardEnum); _assert_(enthalpy_input); - Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input); - Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input); - Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input); - - /*Build friction element, needed later: */ - friction=new Friction(this,3); - - /* Start looping on the number of gauss 2d (nodes on the bedrock) */ - gauss=new GaussPenta(0,1,2,2); - gaussup=new GaussPenta(3,4,5,2); - - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - gaussup->GaussPoint(ig); - - GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); - GetNodalFunctionsP1(&basis[0], gauss); - - enthalpy_input->GetInputValue(&enthalpy,gauss); - pressure_input->GetInputValue(&pressure,gauss); - watercolumn_input->GetInputValue(&watercolumn,gauss); - - if((watercolumn<=0.) && (enthalpyPureIceEnthalpy(pressure))){ - /* the above check is equivalent to - NOT ((watercolumn>0.) AND (enthalpyGetInputValue(&geothermalflux_value,gauss); - - friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vz_input->GetInputValue(&vz,gauss); - basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)+pow(vz,2.0)); - - heatflux_value=(basalfriction+geothermalflux_value)/(rho_ice); - scalar=gauss->weight*Jdet2d*heatflux_value; - if(reCast(dt)) scalar=dt*scalar; - for(i=0;ivalues[i]+=scalar*basis[i]; - } - else if(enthalpy>=matpar->PureIceEnthalpy(pressure)){ - /* check positive thickness of temperate basal ice layer */ - enthalpy_input->GetInputValue(&enthalpyup,gaussup); - pressure_input->GetInputValue(&pressureup,gaussup); - if(enthalpyup>=matpar->PureIceEnthalpy(pressureup)){ - // TODO: temperate ice has positive thickness: grad enthalpy*n=0. - } - else{ - // only base temperate, set Dirichlet BCs in Penta::UpdateBasalConstraintsEnthalpy() - } - } - else { - // base cold, but watercolumn positive. Set base to h_pmp. - } - } - - /*Clean up and return*/ - delete gauss; - delete gaussup; - delete friction; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorMelting {{{*/ -ElementVector* Penta::CreatePVectorMelting(void){ - return NULL; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorThermal {{{*/ -ElementVector* Penta::CreatePVectorThermal(void){ - - /*compute all load vectors for 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 {{{*/ -ElementVector* Penta::CreatePVectorThermalVolume(void){ - - /*Constants*/ - const int numdof=NUMVERTICES*NDOF1; - - /*Intermediaries*/ - int i; - int stabilization; - IssmDouble Jdet,phi,dt; - IssmDouble rho_ice,heatcapacity; - IssmDouble thermalconductivity,kappa; - IssmDouble viscosity,temperature; - IssmDouble tau_parameter,diameter; - IssmDouble u,v,w; - IssmDouble scalar_def,scalar_transient; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble L[numdof]; - IssmDouble dbasis[3][6]; - IssmDouble epsilon[6]; - GaussPenta *gauss=NULL; - - /*Initialize Element vector*/ - ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - rho_ice=matpar->GetRhoIce(); - heatcapacity=matpar->GetHeatCapacity(); - thermalconductivity=matpar->GetThermalConductivity(); - kappa=thermalconductivity/(rho_ice*heatcapacity); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - this->parameters->FindParam(&stabilization,ThermalStabilizationEnum); - 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* temperature_input=NULL; - if(reCast(dt)) temperature_input=inputs->GetInput(TemperatureEnum); _assert_(inputs); - if(stabilization==2) diameter=MinEdgeLength(&xyz_list[0][0]); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(3,3); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctionsP1(&L[0], gauss); - - this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - material->GetViscosity3dFS(&viscosity,&epsilon[0]); - GetPhi(&phi, &epsilon[0], viscosity); - - scalar_def=phi/(rho_ice*heatcapacity)*Jdet*gauss->weight; - if(reCast(dt)) scalar_def=scalar_def*dt; - - for(i=0;ivalues[i]+=scalar_def*L[i]; - - /* Build transient now */ - if(reCast(dt)){ - temperature_input->GetInputValue(&temperature, gauss); - scalar_transient=temperature*Jdet*gauss->weight; - for(i=0;ivalues[i]+=scalar_transient*L[i]; - } - - if(stabilization==2){ - GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss); - - vx_input->GetInputValue(&u, gauss); - vy_input->GetInputValue(&v, gauss); - vz_input->GetInputValue(&w, gauss); - - 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]); - if(reCast(dt)){ - for(i=0;ivalues[i]+=tau_parameter*scalar_transient*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]); - } - } - } - - /*Clean up and return*/ - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorThermalShelf {{{*/ -ElementVector* Penta::CreatePVectorThermalShelf(void){ - - /*Constants*/ - const int numdof=NUMVERTICES*NDOF1; - - /*Intermediaries */ - int i,j; - IssmDouble Jdet2d; - IssmDouble mixed_layer_capacity,thermal_exchange_velocity; - IssmDouble rho_ice,rho_water,pressure,dt,scalar_ocean; - IssmDouble heatcapacity,t_pmp; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_tria[NUMVERTICES2D][3]; - IssmDouble basis[NUMVERTICES]; - GaussPenta* gauss=NULL; - - /* Ice/ocean heat exchange flux on ice shelf base */ - if (!IsOnBed() || !IsFloating()) return NULL; - - /*Initialize Element vector*/ - ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - for(i=0;iGetMixedLayerCapacity(); - thermal_exchange_velocity=matpar->GetThermalExchangeVelocity(); - rho_water=matpar->GetRhoWater(); - rho_ice=matpar->GetRhoIce(); - heatcapacity=matpar->GetHeatCapacity(); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - 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(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); - GetNodalFunctionsP1(&basis[0], gauss); - - pressure_input->GetInputValue(&pressure,gauss); - t_pmp=matpar->TMeltingPoint(pressure); - - scalar_ocean=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(t_pmp)/(heatcapacity*rho_ice); - if(reCast(dt)) scalar_ocean=dt*scalar_ocean; - - for(i=0;ivalues[i]+=scalar_ocean*basis[i]; - } - - /*Clean up and return*/ - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorThermalSheet {{{*/ -ElementVector* Penta::CreatePVectorThermalSheet(void){ - - /*Constants*/ - const int numdof=NUMVERTICES*NDOF1; - - /*Intermediaries */ - int i,j; - int analysis_type; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_tria[NUMVERTICES2D][3]={0.0}; - IssmDouble Jdet2d,dt; - IssmDouble rho_ice,heatcapacity,geothermalflux_value; - IssmDouble basalfriction,alpha2,vx,vy; - IssmDouble basis[NUMVERTICES]; - IssmDouble scalar; - Friction* friction=NULL; - GaussPenta* gauss=NULL; - - /* Geothermal flux on ice sheet base and basal friction */ - if (!IsOnBed() || IsFloating()) return NULL; - - /*Initialize Element vector*/ - ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - for(i=0;iFindParam(&analysis_type,AnalysisTypeEnum); - rho_ice=matpar->GetRhoIce(); - heatcapacity=matpar->GetHeatCapacity(); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - 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(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input); - - /*Build frictoin element, needed later: */ - friction=new Friction(this,3); - - /* Start looping on the number of gauss 2d (nodes on the bedrock) */ - gauss=new GaussPenta(0,1,2,2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); - GetNodalFunctionsP1(&basis[0], gauss); - - geothermalflux_input->GetInputValue(&geothermalflux_value,gauss); - friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)); - - scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice); - if(reCast(dt)) scalar=dt*scalar; - - for(i=0;ivalues[i]+=scalar*basis[i]; - } - - /*Clean up and return*/ - delete gauss; - delete friction; - return pe; -} -/*}}}*/ /*FUNCTION Penta::UpdateBasalConstraintsEnthalpy{{{*/ void Penta::UpdateBasalConstraintsEnthalpy(void){ @@ -5118,427 +4475,4 @@ delete gauss; return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorAdjointHoriz{{{*/ -ElementVector* Penta::CreatePVectorAdjointHoriz(void){ - - int approximation; - inputs->GetInputValue(&approximation,ApproximationEnum); - - switch(approximation){ - case SSAApproximationEnum: - return CreatePVectorAdjointSSA(); - case HOApproximationEnum: - return CreatePVectorAdjointHO(); - case NoneApproximationEnum: - return NULL; - case FSApproximationEnum: - return CreatePVectorAdjointFS(); - default: - _error_("Approximation " << EnumToStringx(approximation) << " not supported yet"); - } -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorAdjointSSA{{{*/ -ElementVector* Penta::CreatePVectorAdjointSSA(){ - - if (!IsOnBed()) return NULL; - - /*Call Tria function*/ - Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. - ElementVector* pe=tria->CreatePVectorAdjointHoriz(); - delete tria->material; delete tria; - - /*clean up and return*/ - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorAdjointHO{{{*/ -ElementVector* Penta::CreatePVectorAdjointHO(void){ - - /*Nothing to be done if not on surface*/ - if (!IsOnSurface()) return NULL; - - /*Intermediaries */ - int i,j,resp; - int *responses=NULL; - int num_responses; - IssmDouble Jdet2d; - IssmDouble obs_velocity_mag,velocity_mag; - IssmDouble dux,duy; - IssmDouble epsvel=2.220446049250313e-16; - IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/ - IssmDouble scalex=0.,scaley=0.,scale=0.,S=0.; - IssmDouble vx,vy,vxobs,vyobs,weight; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_tria[NUMVERTICES2D][3]; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element vector*/ - ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum); - this->parameters->FindParam(&responses,NULL,InversionCostFunctionsEnum); - Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); - Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input =inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vxobs_input =inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input); - Input* vyobs_input =inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input); - - for(i=0;iGetInputValue(&S,SurfaceAreaEnum); break; - } - } - - /* Start looping on the number of gaussian points: */ - GaussPenta* gauss=new GaussPenta(3,4,5,4); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - /* Get Jacobian determinant: */ - GetTriaJacobianDeterminant(&Jdet2d,&xyz_list_tria[0][0],gauss); - - /*Get all parameters at gaussian point*/ - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vxobs_input->GetInputValue(&vxobs,gauss); - vyobs_input->GetInputValue(&vyobs,gauss); - GetNodalFunctions(basis, gauss); - - /*Loop over all requested responses*/ - for(resp=0;respGetInputValue(&weight,gauss,responses[resp]); - - switch(responses[resp]){ - case SurfaceAbsVelMisfitEnum: - /* - * 1 [ 2 2 ] - * J = --- | (u - u ) + (v - v ) | - * 2 [ obs obs ] - * - * dJ - * DU = - -- = (u - u ) - * du obs - */ - for(i=0;ivalues[i*NDOF2+0]+=dux*weight*Jdet2d*gauss->weight*basis[i]; - pe->values[i*NDOF2+1]+=duy*weight*Jdet2d*gauss->weight*basis[i]; - } - break; - case SurfaceRelVelMisfitEnum: - /* - * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] - * J = --- | ------------- (u - u ) + ------------- (v - v ) | - * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] - * obs obs - * - * dJ \bar{v}^2 - * DU = - -- = ------------- (u - u ) - * du (u + eps)^2 obs - * obs - */ - for(i=0;ivalues[i*NDOF2+0]+=dux*weight*Jdet2d*gauss->weight*basis[i]; - pe->values[i*NDOF2+1]+=duy*weight*Jdet2d*gauss->weight*basis[i]; - } - break; - case SurfaceLogVelMisfitEnum: - /* - * [ vel + eps ] 2 - * J = 4 \bar{v}^2 | log ( ----------- ) | - * [ vel + eps ] - * obs - * - * dJ 2 * log(...) - * DU = - -- = - 4 \bar{v}^2 ------------- u - * du vel^2 + eps - * - */ - for(i=0;ivalues[i*NDOF2+0]+=dux*weight*Jdet2d*gauss->weight*basis[i]; - pe->values[i*NDOF2+1]+=duy*weight*Jdet2d*gauss->weight*basis[i]; - } - break; - case SurfaceAverageVelMisfitEnum: - /* - * 1 2 2 - * J = --- sqrt( (u - u ) + (v - v ) ) - * S obs obs - * - * dJ 1 1 - * DU = - -- = - --- ----------- * 2 (u - u ) - * du S 2 sqrt(...) obs - */ - for(i=0;ivalues[i*NDOF2+0]+=dux*weight*Jdet2d*gauss->weight*basis[i]; - pe->values[i*NDOF2+1]+=duy*weight*Jdet2d*gauss->weight*basis[i]; - } - break; - case SurfaceLogVxVyMisfitEnum: - /* - * 1 [ |u| + eps 2 |v| + eps 2 ] - * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | - * 2 [ |u |+ eps |v |+ eps ] - * obs obs - * dJ 1 u 1 - * DU = - -- = - \bar{v}^2 log(u...) --------- ---- ~ - \bar{v}^2 log(u...) ------ - * du |u| + eps |u| u + eps - */ - for(i=0;ivalues[i*NDOF2+0]+=dux*weight*Jdet2d*gauss->weight*basis[i]; - pe->values[i*NDOF2+1]+=duy*weight*Jdet2d*gauss->weight*basis[i]; - } - break; - case DragCoefficientAbsGradientEnum: - /*Nothing in P vector*/ - break; - case ThicknessAbsGradientEnum: - /*Nothing in P vector*/ - break; - case ThicknessAlongGradientEnum: - /*Nothing in P vector*/ - break; - case ThicknessAcrossGradientEnum: - /*Nothing in P vector*/ - break; - case RheologyBbarAbsGradientEnum: - /*Nothing in P vector*/ - break; - default: - _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); - } - } - } - - /*Clean up and return*/ - xDelete(basis); - xDelete(responses); - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorAdjointFS{{{*/ -ElementVector* Penta::CreatePVectorAdjointFS(void){ - - if(!IsOnSurface()) return NULL; - - /*Intermediaries */ - int i,j,resp; - int *responses= NULL; - int num_responses; - IssmDouble Jdet2d; - IssmDouble obs_velocity_mag,velocity_mag; - IssmDouble dux,duy; - IssmDouble epsvel = 2.220446049250313e-16; - IssmDouble meanvel = 3.170979198376458e-05; /*1000 m/yr */ - IssmDouble scalex = 0,scaley=0,scale=0,S=0; - IssmDouble vx,vy,vxobs,vyobs,weight; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_tria[NUMVERTICES2D][3]; - - /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = this->NumberofNodesVelocity(); - int pnumnodes = this->NumberofNodesPressure(); - - /*Prepare coordinate system list*/ - int* cs_list = xNew(vnumnodes+pnumnodes); - for(i=0;iparameters,FSApproximationEnum); - IssmDouble* vbasis = xNew(vnumnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum); - this->parameters->FindParam(&responses,NULL,InversionCostFunctionsEnum); - Input* weights_input = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); - Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vxobs_input = inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input); - Input* vyobs_input = inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input); - - for(i=0;iGetInputValue(&S,SurfaceAreaEnum); break; - } - } - - /* Start looping on the number of gaussian points: */ - GaussPenta* gauss=new GaussPenta(3,4,5,4); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - /* Get Jacobian determinant: */ - GetTriaJacobianDeterminant(&Jdet2d,&xyz_list_tria[0][0],gauss); - - /*Get all parameters at gaussian point*/ - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vxobs_input->GetInputValue(&vxobs,gauss); - vyobs_input->GetInputValue(&vyobs,gauss); - GetNodalFunctionsVelocity(vbasis,gauss); - - /*Loop over all requested responses*/ - for(resp=0;respGetInputValue(&weight,gauss,responses[resp]); - - switch(responses[resp]){ - - case SurfaceAbsVelMisfitEnum: - /* - * 1 [ 2 2 ] - * J = --- | (u - u ) + (v - v ) | - * 2 [ obs obs ] - * - * dJ - * DU = - -- = (u - u ) - * du obs - */ - for(i=0;ivalues[i*NDOF3+0]+=dux*weight*Jdet2d*gauss->weight*vbasis[i]; - pe->values[i*NDOF3+1]+=duy*weight*Jdet2d*gauss->weight*vbasis[i]; - } - break; - case SurfaceRelVelMisfitEnum: - /* - * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] - * J = --- | ------------- (u - u ) + ------------- (v - v ) | - * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] - * obs obs - * - * dJ \bar{v}^2 - * DU = - -- = ------------- (u - u ) - * du (u + eps)^2 obs - * obs - */ - for(i=0;ivalues[i*NDOF3+0]+=dux*weight*Jdet2d*gauss->weight*vbasis[i]; - pe->values[i*NDOF3+1]+=duy*weight*Jdet2d*gauss->weight*vbasis[i]; - } - break; - case SurfaceLogVelMisfitEnum: - /* - * [ vel + eps ] 2 - * J = 4 \bar{v}^2 | log ( ----------- ) | - * [ vel + eps ] - * obs - * - * dJ 2 * log(...) - * DU = - -- = - 4 \bar{v}^2 ------------- u - * du vel^2 + eps - * - */ - for(i=0;ivalues[i*NDOF3+0]+=dux*weight*Jdet2d*gauss->weight*vbasis[i]; - pe->values[i*NDOF3+1]+=duy*weight*Jdet2d*gauss->weight*vbasis[i]; - } - break; - case SurfaceAverageVelMisfitEnum: - /* - * 1 2 2 - * J = --- sqrt( (u - u ) + (v - v ) ) - * S obs obs - * - * dJ 1 1 - * DU = - -- = - --- ----------- * 2 (u - u ) - * du S 2 sqrt(...) obs - */ - for(i=0;ivalues[i*NDOF3+0]+=dux*weight*Jdet2d*gauss->weight*vbasis[i]; - pe->values[i*NDOF3+1]+=duy*weight*Jdet2d*gauss->weight*vbasis[i]; - } - break; - case SurfaceLogVxVyMisfitEnum: - /* - * 1 [ |u| + eps 2 |v| + eps 2 ] - * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | - * 2 [ |u |+ eps |v |+ eps ] - * obs obs - * dJ 1 u 1 - * DU = - -- = - \bar{v}^2 log(u...) --------- ---- ~ - \bar{v}^2 log(u...) ------ - * du |u| + eps |u| u + eps - */ - for(i=0;ivalues[i*NDOF3+0]+=dux*weight*Jdet2d*gauss->weight*vbasis[i]; - pe->values[i*NDOF3+1]+=duy*weight*Jdet2d*gauss->weight*vbasis[i]; - } - break; - case DragCoefficientAbsGradientEnum: - /*Nothing in P vector*/ - break; - case ThicknessAbsGradientEnum: - /*Nothing in P vector*/ - break; - case ThicknessAcrossGradientEnum: - /*Nothing in P vector*/ - break; - case ThicknessAlongGradientEnum: - /*Nothing in P vector*/ - break; - case RheologyBbarAbsGradientEnum: - /*Nothing in P vector*/ - break; - default: - _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); - } - } - } - - /*Clean up and return*/ - xDelete(responses); - xDelete(cs_list); - xDelete(vbasis); - delete gauss; - return pe; } /*}}}*/ @@ -7867,1150 +6801,4 @@ } /*}}}*/ -/*FUNCTION Penta::CreatePVectorCouplingSSAFS {{{*/ -ElementVector* Penta::CreatePVectorCouplingSSAFS(void){ - - /*compute all load vectors for this element*/ - ElementVector* pe1=CreatePVectorCouplingSSAFSViscous(); - ElementVector* pe2=CreatePVectorCouplingSSAFSFriction(); - ElementVector* pe =new ElementVector(pe1,pe2); - - /*clean-up and return*/ - delete pe1; - delete pe2; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorCouplingSSAFSViscous {{{*/ -ElementVector* Penta::CreatePVectorCouplingSSAFSViscous(void){ - - /*Intermediaries */ - int i,approximation; - IssmDouble viscosity,Jdet; - IssmDouble FSreconditioning; - IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - IssmDouble dw[3]; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble basis[6]; //for the six nodes of the penta - IssmDouble dbasis[3][6]; //for the six nodes of the penta - GaussPenta *gauss=NULL; - - /*Initialize Element vector and return if necessary*/ - inputs->GetInputValue(&approximation,ApproximationEnum); - if(approximation!=SSAFSApproximationEnum) return NULL; - int vnumnodes = this->NumberofNodesVelocity(); - int pnumnodes = this->NumberofNodesPressure(); - - /*Prepare coordinate system list*/ - int* cs_list = xNew(vnumnodes+pnumnodes); - for(i=0;iparameters,FSvelocityEnum); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); - 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* vzSSA_input=inputs->GetInput(VzSSAEnum); _assert_(vzSSA_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(5,5); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctionsP1(&basis[0], gauss); - GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss); - - vzSSA_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss); - - this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - material->GetViscosity3dFS(&viscosity,&epsilon[0]); - - for(i=0;ivalues[i*NDOF3+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i]; - pe->values[i*NDOF3+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i]; - pe->values[i*NDOF3+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]); - pe->values[NDOF3*vnumnodes+i]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i]; - } - } - - /*Transform coordinate system*/ - ::TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list); - - /*Clean up and return*/ - xDelete(cs_list); - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorCouplingSSAFSFriction{{{*/ -ElementVector* Penta::CreatePVectorCouplingSSAFSFriction(void){ - - /*Intermediaries*/ - int i,j; - int approximation,analysis_type; - IssmDouble Jdet,Jdet2d; - IssmDouble FSreconditioning; - IssmDouble bed_normal[3]; - IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - IssmDouble viscosity, w, alpha2_gauss; - IssmDouble dw[3]; - IssmDouble xyz_list_tria[NUMVERTICES2D][3]; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble basis[6]; //for the six nodes of the penta - Friction* friction=NULL; - GaussPenta *gauss=NULL; - - /*Initialize Element vector and return if necessary*/ - if(!IsOnBed() || IsFloating()) return NULL; - inputs->GetInputValue(&approximation,ApproximationEnum); - if(approximation!=SSAFSApproximationEnum) return NULL; - int vnumnodes = this->NumberofNodesVelocity(); - int pnumnodes = this->NumberofNodesPressure(); - - /*Prepare coordinate system list*/ - int* cs_list = xNew(vnumnodes+pnumnodes); - for(i=0;iparameters,FSvelocityEnum); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); - 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* vzSSA_input=inputs->GetInput(VzSSAEnum); _assert_(vzSSA_input); - - for(i=0;ibegin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); - GetNodalFunctionsP1(basis, gauss); - - vzSSA_input->GetInputValue(&w, gauss); - vzSSA_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss); - - NormalBase(&bed_normal[0],&xyz_list_tria[0][0]); - this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - material->GetViscosity3dFS(&viscosity,&epsilon[0]); - friction->GetAlpha2(&alpha2_gauss, gauss,vx_input,vy_input,vz_input); - - for(i=0;ivalues[i*NDOF3+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i]; - pe->values[i*NDOF3+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i]; - pe->values[i*NDOF3+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*basis[i]; - } - } - - /*Transform coordinate system*/ - ::TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list); - - /*Clean up and return*/ - xDelete(cs_list); - delete gauss; - delete friction; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorCouplingHOFS {{{*/ -ElementVector* Penta::CreatePVectorCouplingHOFS(void){ - - /*compute all load vectors for this element*/ - ElementVector* pe1=CreatePVectorCouplingHOFSViscous(); - ElementVector* pe2=CreatePVectorCouplingHOFSFriction(); - ElementVector* pe =new ElementVector(pe1,pe2); - - /*clean-up and return*/ - delete pe1; - delete pe2; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorCouplingHOFSViscous {{{*/ -ElementVector* Penta::CreatePVectorCouplingHOFSViscous(void){ - - /*Intermediaries */ - int i; - int approximation; - IssmDouble viscosity,Jdet; - IssmDouble FSreconditioning; - IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - IssmDouble dw[3]; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble basis[6]; //for the six nodes of the penta - IssmDouble dbasis[3][6]; //for the six nodes of the penta - GaussPenta *gauss=NULL; - - /*Initialize Element vector and return if necessary*/ - inputs->GetInputValue(&approximation,ApproximationEnum); - if(approximation!=HOFSApproximationEnum) return NULL; - int vnumnodes = this->NumberofNodesVelocity(); - int pnumnodes = this->NumberofNodesPressure(); - - /*Prepare coordinate system list*/ - int* cs_list = xNew(vnumnodes+pnumnodes); - for(i=0;iparameters,FSvelocityEnum); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); - 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* vzHO_input=inputs->GetInput(VzHOEnum); _assert_(vzHO_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(5,5); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctionsP1(&basis[0], gauss); - GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss); - - vzHO_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss); - - this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - material->GetViscosity3dFS(&viscosity,&epsilon[0]); - - for(i=0;ivalues[i*NDOF3+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i]; - pe->values[i*NDOF3+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i]; - pe->values[i*NDOF3+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]); - pe->values[NDOF3*vnumnodes+i]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i]; - } - } - - /*Transform coordinate system*/ - ::TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list); - - /*Clean up and return*/ - xDelete(cs_list); - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorCouplingHOFSFriction{{{*/ -ElementVector* Penta::CreatePVectorCouplingHOFSFriction(void){ - - /*Intermediaries*/ - int i,j; - int approximation,analysis_type; - IssmDouble Jdet,Jdet2d; - IssmDouble FSreconditioning; - IssmDouble bed_normal[3]; - IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - IssmDouble viscosity, w, alpha2_gauss; - IssmDouble dw[3]; - IssmDouble xyz_list_tria[NUMVERTICES2D][3]; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble basis[6]; //for the six nodes of the penta - Friction* friction=NULL; - GaussPenta *gauss=NULL; - - /*Initialize Element vector and return if necessary*/ - if(!IsOnBed() || IsFloating()) return NULL; - inputs->GetInputValue(&approximation,ApproximationEnum); - if(approximation!=HOFSApproximationEnum) return NULL; - - int vnumnodes = this->NumberofNodesVelocity(); - int pnumnodes = this->NumberofNodesPressure(); - int numnodes = vnumnodes+pnumnodes; - - /*Prepare coordinate system list*/ - int* cs_list = xNew(vnumnodes+pnumnodes); - for(i=0;iparameters,FSvelocityEnum); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); - 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* vzHO_input=inputs->GetInput(VzHOEnum); _assert_(vzHO_input); - - for(i=0;ibegin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); - GetNodalFunctionsP1(basis, gauss); - - vzHO_input->GetInputValue(&w, gauss); - vzHO_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss); - - NormalBase(&bed_normal[0],&xyz_list_tria[0][0]); - this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - material->GetViscosity3dFS(&viscosity,&epsilon[0]); - friction->GetAlpha2(&alpha2_gauss, gauss,vx_input,vy_input,vz_input); - - for(i=0;ivalues[i*NDOF3+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i]; - pe->values[i*NDOF3+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i]; - pe->values[i*NDOF3+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*basis[i]; - } - } - - /*Transform coordinate system*/ - ::TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list); - - /*Clean up and return*/ - xDelete(cs_list); - delete gauss; - delete friction; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorStressbalanceHoriz{{{*/ -ElementVector* Penta::CreatePVectorStressbalanceHoriz(void){ - - int approximation; - inputs->GetInputValue(&approximation,ApproximationEnum); - - switch(approximation){ - case SSAApproximationEnum: - return CreatePVectorStressbalanceSSA(); - case HOApproximationEnum: - return CreatePVectorStressbalanceHO(); - case L1L2ApproximationEnum: - return CreatePVectorStressbalanceL1L2(); - case SIAApproximationEnum: - return NULL; - case NoneApproximationEnum: - return NULL; - case FSApproximationEnum: - return CreatePVectorStressbalanceFS(); - case SSAHOApproximationEnum: - return CreatePVectorStressbalanceSSAHO(); - case SSAFSApproximationEnum: - return CreatePVectorStressbalanceSSAFS(); - case HOFSApproximationEnum: - return CreatePVectorStressbalanceHOFS(); - default: - _error_("Approximation " << EnumToStringx(approximation) << " not supported yet"); - } -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorStressbalanceSSAHO{{{*/ -ElementVector* Penta::CreatePVectorStressbalanceSSAHO(void){ - - /*compute all load vectors for this element*/ - ElementVector* pe1=CreatePVectorStressbalanceSSA(); - ElementVector* pe2=CreatePVectorStressbalanceHO(); - ElementVector* pe =new ElementVector(pe1,pe2); - - /*clean-up and return*/ - delete pe1; - delete pe2; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorStressbalanceSSAFS{{{*/ -ElementVector* Penta::CreatePVectorStressbalanceSSAFS(void){ - - /*compute all load vectors for this element*/ - int init = this->element_type; - this->element_type=P1Enum; - ElementVector* pe1=CreatePVectorStressbalanceSSA(); - this->element_type=init; - ElementVector* pe2=CreatePVectorStressbalanceFS(); - int indices[3]={18,19,20}; - this->element_type=MINIcondensedEnum; - ElementMatrix* Ke = CreateKMatrixStressbalanceFS(); - this->element_type=init; - pe2->StaticCondensation(Ke,3,&indices[0]); - delete Ke; - ElementVector* pe3=CreatePVectorCouplingSSAFS(); - ElementVector* pe =new ElementVector(pe1,pe2,pe3); - - /*clean-up and return*/ - delete pe1; - delete pe2; - delete pe3; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorStressbalanceHOFS{{{*/ -ElementVector* Penta::CreatePVectorStressbalanceHOFS(void){ - - /*compute all load vectors for this element*/ - int init = this->element_type; - this->element_type=P1Enum; - ElementVector* pe1=CreatePVectorStressbalanceHO(); - this->element_type=init; - ElementVector* pe2=CreatePVectorStressbalanceFS(); - int indices[3]={18,19,20}; - this->element_type=MINIcondensedEnum; - ElementMatrix* Ke = CreateKMatrixStressbalanceFS(); - this->element_type=init; - pe2->StaticCondensation(Ke,3,&indices[0]); - delete Ke; - ElementVector* pe3=CreatePVectorCouplingHOFS(); - ElementVector* pe =new ElementVector(pe1,pe2,pe3); - - /*clean-up and return*/ - delete pe1; - delete pe2; - delete pe3; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorStressbalanceSIA{{{*/ -ElementVector* Penta::CreatePVectorStressbalanceSIA(void){ - - /*Intermediaries*/ - int i,j; - int node0,node1; - IssmDouble connectivity[2]; - IssmDouble Jdet; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_segment[2][3]; - IssmDouble z_list[NUMVERTICES]; - IssmDouble slope[2]; - IssmDouble slope2,constant_part,drag; - IssmDouble rho_ice,gravity,n,B; - IssmDouble ub,vb,z_g,surface,thickness; - GaussPenta* gauss=NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); _assert_(numnodes==6); - - /*Initialize Element vector*/ - ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - rho_ice=matpar->GetRhoIce(); - gravity=matpar->GetG(); - n=material->GetN(); - B=material->GetB(); - Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); - Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input); - Input* slopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(slopex_input); - Input* slopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(slopey_input); - Input* drag_input=inputs->GetInput(FrictionCoefficientEnum); _assert_(drag_input); - for(i=0;iConnectivity(); - connectivity[1]=(IssmDouble)vertices[node1]->Connectivity(); - - /*Loop on the Gauss points: */ - gauss=new GaussPenta(node0,node1,3); - for(int ig=gauss->begin();igend();ig++){ - gauss->GaussPoint(ig); - - slopex_input->GetInputValue(&slope[0],gauss); - slopey_input->GetInputValue(&slope[1],gauss); - surface_input->GetInputValue(&surface,gauss); - thickness_input->GetInputValue(&thickness,gauss); - - slope2=pow(slope[0],2)+pow(slope[1],2); - constant_part=-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.)); - - PentaRef::GetInputValue(&z_g,&z_list[0],gauss); - GetSegmentJacobianDeterminant(&Jdet,&xyz_list_segment[0][0],gauss); - - if(IsOnSurface()){ - pe->values[2*node1+0]+=constant_part*pow((surface-z_g)/B,n)*slope[0]*Jdet*gauss->weight/connectivity[1]; - pe->values[2*node1+1]+=constant_part*pow((surface-z_g)/B,n)*slope[1]*Jdet*gauss->weight/connectivity[1]; - } - else{/*connectivity is too large, should take only half on it*/ - pe->values[2*node1+0]+=constant_part*pow((surface-z_g)/B,n)*slope[0]*Jdet*gauss->weight*2./connectivity[1]; - pe->values[2*node1+1]+=constant_part*pow((surface-z_g)/B,n)*slope[1]*Jdet*gauss->weight*2./connectivity[1]; - } - } - - /*Deal with lower surface*/ - if(IsOnBed()){ - drag_input->GetInputValue(&drag,gauss); - - /*Payne 1995 (m = 1, B = 5e-3/yts = 1.58e-10 m/s/Pa*/ - ub=-1.58*1.e-10*rho_ice*gravity*thickness*slope[0]; - vb=-1.58*1.e-10*rho_ice*gravity*thickness*slope[1]; - ///*Ritz et al. 1996*/ - //ub=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[0]/sqrt(slope2); - //vb=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[1]/sqrt(slope2); - ///*Rutt et al. 2009*/ - //ub=-drag*rho_ice*gravity*thickness*slope[0]; - //vb=-drag*rho_ice*gravity*thickness*slope[1]; - - pe->values[2*node0+0]+=ub/connectivity[0]; - pe->values[2*node0+1]+=vb/connectivity[0]; - } - delete gauss; - } - - /*Clean up and return*/ - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorStressbalanceSSA{{{*/ -ElementVector* Penta::CreatePVectorStressbalanceSSA(void){ - - if (!IsOnBed()) return NULL; - - /*Call Tria function*/ - Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. - ElementVector* pe=tria->CreatePVectorStressbalanceSSA(); - delete tria->material; delete tria; - - /*Clean up and return*/ - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorStressbalanceL1L2{{{*/ -ElementVector* Penta::CreatePVectorStressbalanceL1L2(void){ - - if (!IsOnBed()) return NULL; - - /*Call Tria function*/ - Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. - ElementVector* pe=tria->CreatePVectorStressbalanceSSA(); - delete tria->material; delete tria; - - /*Clean up and return*/ - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorStressbalanceHO{{{*/ -ElementVector* Penta::CreatePVectorStressbalanceHO(void){ - - /*compute all load vectors for this element*/ - ElementVector* pe1=CreatePVectorStressbalanceHODrivingStress(); - ElementVector* pe2=CreatePVectorStressbalanceHOFront(); - ElementVector* pe =new ElementVector(pe1,pe2); - - /*clean-up and return*/ - delete pe1; - delete pe2; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorStressbalanceHODrivingStress{{{*/ -ElementVector* Penta::CreatePVectorStressbalanceHODrivingStress(void){ - - /*Intermediaries*/ - IssmDouble Jdet; - IssmDouble slope[3]; //do not put 2! this goes into GetInputDerivativeValue, which addresses slope[3] also! - IssmDouble driving_stress_baseline,thickness; - IssmDouble xyz_list[NUMVERTICES][3]; - GaussPenta *gauss=NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize vector*/ - ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters,HOApproximationEnum); - IssmDouble* basis = xNewZeroInit(numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); - Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(3,3); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - - thickness_input->GetInputValue(&thickness, gauss); - surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss); - - driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG(); - - for(int i=0;ivalues[i*NDOF2+0]+= -driving_stress_baseline*slope[0]*Jdet*gauss->weight*basis[i]; - pe->values[i*NDOF2+1]+= -driving_stress_baseline*slope[1]*Jdet*gauss->weight*basis[i]; - } - } - - /*Transform coordinate system*/ - ::TransformLoadVectorCoord(pe,nodes,numnodes,XYEnum); - - /*Clean up and return*/ - xDelete(basis); - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorStressbalanceHOFront{{{*/ -ElementVector* Penta::CreatePVectorStressbalanceHOFront(void){ - - /*If no front, return NULL*/ - if(!IsZeroLevelset(MaskIceLevelsetEnum)) return NULL; - - /*Intermediaries*/ - IssmDouble rho_ice,rho_water,gravity; - IssmDouble Jdet,surface,z_g,water_pressure,ice_pressure,air_pressure; - IssmDouble surface_under_water,base_under_water,pressure; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble* xyz_list_front = NULL; - IssmDouble area_coordinates[4][3]; - IssmDouble normal[3]; - GaussPenta* gauss; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element vector and other vectors*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters,HOApproximationEnum); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input); - rho_water=matpar->GetRhoWater(); - rho_ice =matpar->GetRhoIce(); - gravity =matpar->GetG(); - ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum); - GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],4); - NormalSection(&normal[0],xyz_list_front); - - /*Initialize gauss points*/ - IssmDouble zmax=xyz_list[0][2]; for(int i=1;i<6;i++) if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2]; - IssmDouble zmin=xyz_list[0][2]; for(int i=1;i<6;i++) if(xyz_list[i][2]0 && zmin<0) gauss=new GaussPenta(area_coordinates,3,10); //refined in vertical because of the sea level discontinuity - else gauss=new GaussPenta(area_coordinates,3,3); - - /* Start looping on the number of gaussian points: */ - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - surface_input->GetInputValue(&surface,gauss); - z_g=GetZcoord(gauss); - GetNodalFunctions(basis,gauss); - GetQuadJacobianDeterminant(&Jdet,xyz_list_front,gauss); - - water_pressure = rho_water*gravity*min(0.,z_g);//0 if the gaussian point is above water level - ice_pressure = rho_ice*gravity*(surface-z_g); - air_pressure = 0; - pressure = ice_pressure + water_pressure + air_pressure; - - for (int i=0;ivalues[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i]; - pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i]; - } - } - - /*Transform coordinate system*/ - ::TransformLoadVectorCoord(pe,nodes,numnodes,XYEnum); - - /*Clean up and return*/ - xDelete(basis); - xDelete(xyz_list_front); - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorStressbalanceFS {{{*/ -ElementVector* Penta::CreatePVectorStressbalanceFS(void){ - - ElementVector* pe1; - ElementVector* pe2; - ElementVector* pe3; - ElementVector* pe; - - /*compute all stiffness matrices for this element*/ - pe1=CreatePVectorStressbalanceFSViscous(); - pe2=CreatePVectorStressbalanceFSShelf(); - pe3=CreatePVectorStressbalanceFSFront(); - pe =new ElementVector(pe1,pe2,pe3); - - /*clean-up and return*/ - delete pe1; - delete pe2; - delete pe3; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorStressbalanceFSFront{{{*/ -ElementVector* Penta::CreatePVectorStressbalanceFSFront(void){ - - /*If no front, return NULL*/ - if(!IsZeroLevelset(MaskIceLevelsetEnum)) return NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int i; - IssmDouble rho_ice,rho_water,gravity; - IssmDouble Jdet,z_g,water_pressure,air_pressure; - IssmDouble surface_under_water,base_under_water,pressure; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble* xyz_list_front = NULL; - IssmDouble area_coordinates[4][3]; - IssmDouble normal[3]; - GaussPenta* gauss; - - /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = this->NumberofNodesVelocity(); - int pnumnodes = this->NumberofNodesPressure(); - - /*Prepare coordinate system list*/ - int* cs_list = xNew(vnumnodes+pnumnodes); - for(i=0;iparameters,FSvelocityEnum); - IssmDouble* vbasis = xNew(vnumnodes); - - /*Retrieve all inputs and parameters*/ - rho_water=matpar->GetRhoWater(); - rho_ice =matpar->GetRhoIce(); - gravity =matpar->GetG(); - ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum); - GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],4); - NormalSection(&normal[0],xyz_list_front); - - /*Start looping on Gaussian points*/ - IssmDouble zmax=xyz_list[0][2]; for(int i=1;i<6;i++) if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2]; - IssmDouble zmin=xyz_list[0][2]; for(int i=1;i<6;i++) if(xyz_list[i][2]0 && zmin<0) gauss=new GaussPenta(area_coordinates,3,30); //refined in vertical because of the sea level discontinuity - else gauss=new GaussPenta(area_coordinates,3,3); - - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - z_g=GetZcoord(gauss); - GetNodalFunctionsVelocity(vbasis,gauss); - GetQuadJacobianDeterminant(&Jdet,xyz_list_front,gauss); - - water_pressure=rho_water*gravity*min(0.,z_g);//0 if the gaussian point is above water level - air_pressure=0; - pressure = water_pressure + air_pressure; - - for(i=0;ivalues[3*i+0]+= pressure*Jdet*gauss->weight*normal[0]*vbasis[i]; - pe->values[3*i+1]+= pressure*Jdet*gauss->weight*normal[1]*vbasis[i]; - pe->values[3*i+2]+= pressure*Jdet*gauss->weight*normal[2]*vbasis[i]; - } - } - - /*Transform coordinate system*/ - ::TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list); - - /*Clean up and return*/ - xDelete(cs_list); - xDelete(xyz_list_front); - xDelete(vbasis); - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::PVectorGLSstabilization{{{*/ -void Penta::PVectorGLSstabilization(ElementVector* pe){ - - /*Intermediaries*/ - int i,j; - IssmDouble Jdet,gravity,rho_ice,B,D_scalar_stab,viscosity; - IssmDouble forcex,forcey,forcez,diameter,FSreconditioning; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - GaussPenta *gauss=NULL; - - /*Stabilization*/ - IssmDouble dbasis[3][6]; - IssmDouble dmu[3]; - IssmDouble mu; - IssmDouble dnodalbasis[6][6][3]; - IssmDouble SW[6][4][4]; - int p,ii; - int c=3; //index of pressure - - /*Retrieve all inputs and parameters*/ - rho_ice=matpar->GetRhoIce(); - gravity=matpar->GetG(); - B=material->GetB(); - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - 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); - - /*Find minimal length*/ - diameter=MinEdgeLength(&xyz_list[0][0]); - - gauss=new GaussPenta(); - for(int iv=0;iv<6;iv++){ - gauss->GaussVertex(iv); - GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss); - for(i=0;i<6;i++){ - for(j=0;j<3;j++){ - dnodalbasis[i][iv][j] = dbasis[j][i]; - } - } - } - delete gauss; - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(5,5); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - material->GetViscosity3dFS(&viscosity,&epsilon[0]); - - GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss); - dmu[0]=0.; dmu[1]=0.; dmu[2]=0.; - mu = 2.*viscosity; - for(p=0;p<6;p++) for(i=0;i<4;i++) for(j=0;j<4;j++){ - SW[p][i][j]=0.; - } - for(p=0;p<6;p++){ - for(i=0;i<3;i++){ - SW[p][c][i] += dbasis[i][p]; - for(j=0;j<3;j++){ - SW[p][i][i] += -dmu[j]*dbasis[j][p]; - SW[p][j][i] += -dmu[j]*dbasis[i][p]; - for(ii=0;ii<6;ii++){ - SW[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii]; - SW[p][j][i] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii]; - } - } - } - } - IssmDouble tau = 1./3.*pow(diameter,2)/(8.*2.*viscosity); - for(p=0;p<6;p++){ - for(j=0;j<4;j++){ - pe->values[p*4+j] += gauss->weight*Jdet*tau*rho_ice*forcex*SW[p][j][0]; - pe->values[p*4+j] += gauss->weight*Jdet*tau*rho_ice*forcey*SW[p][j][1]; - pe->values[p*4+j] += gauss->weight*Jdet*tau*rho_ice*(forcez-gravity)*SW[p][j][2]; - } - } - } - - /*Clean up*/ - delete gauss; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorStressbalanceFSViscous {{{*/ -ElementVector* Penta::CreatePVectorStressbalanceFSViscous(void){ - - /*Intermediaries*/ - int i; - int approximation; - IssmDouble Jdet,gravity,rho_ice; - IssmDouble forcex,forcey,forcez; - IssmDouble xyz_list[NUMVERTICES][3]; - GaussPenta *gauss=NULL; - - /*Initialize Element vector and return if necessary*/ - inputs->GetInputValue(&approximation,ApproximationEnum); - if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = this->NumberofNodesVelocity(); - int pnumnodes = this->NumberofNodesPressure(); - - /*Prepare coordinate system list*/ - int* cs_list = xNew(vnumnodes+pnumnodes); - for(i=0;iparameters,FSvelocityEnum); - IssmDouble* vbasis = xNew(vnumnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* loadingforcex_input=inputs->GetInput(LoadingforceXEnum); _assert_(loadingforcex_input); - Input* loadingforcey_input=inputs->GetInput(LoadingforceYEnum); _assert_(loadingforcey_input); - Input* loadingforcez_input=inputs->GetInput(LoadingforceZEnum); _assert_(loadingforcez_input); - rho_ice = matpar->GetRhoIce(); - gravity = matpar->GetG(); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(5,5); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctionsVelocity(vbasis, gauss); - - loadingforcex_input->GetInputValue(&forcex,gauss); - loadingforcey_input->GetInputValue(&forcey,gauss); - loadingforcez_input->GetInputValue(&forcez,gauss); - - for(i=0;ivalues[i*NDOF3+2] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i]; - pe->values[i*NDOF3+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i]; - pe->values[i*NDOF3+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i]; - pe->values[i*NDOF3+2] += +rho_ice*forcez *Jdet*gauss->weight*vbasis[i]; - } - } - - /*Transform coordinate system*/ - ::TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list); - - /*Clean up and return*/ - delete gauss; - xDelete(cs_list); - xDelete(vbasis); - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorStressbalanceFSShelf{{{*/ -ElementVector* Penta::CreatePVectorStressbalanceFSShelf(void){ - - /*Intermediaries*/ - int i,j; - int approximation,shelf_dampening; - IssmDouble gravity,rho_water,bed,water_pressure; - IssmDouble damper,normal_vel,vx,vy,vz,dt; - IssmDouble xyz_list_tria[NUMVERTICES2D][3]; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble bed_normal[3]; - IssmDouble dz[3]; - IssmDouble Jdet2d; - GaussPenta *gauss=NULL; - - /*Initialize Element vector and return if necessary*/ - if(!IsOnBed() || !IsFloating()) return NULL; - inputs->GetInputValue(&approximation,ApproximationEnum); - if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = this->NumberofNodesVelocity(); - int pnumnodes = this->NumberofNodesPressure(); - - /*Prepare coordinate system list*/ - int* cs_list = xNew(vnumnodes+pnumnodes); - for(i=0;iparameters,FSvelocityEnum); - IssmDouble* vbasis = xNew(vnumnodes); - - /*Retrieve all inputs and parameters*/ - this->parameters->FindParam(&shelf_dampening,StressbalanceShelfDampeningEnum); - rho_water=matpar->GetRhoWater(); - gravity=matpar->GetG(); - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* bed_input=inputs->GetInput(BedEnum); _assert_(bed_input); - 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); - - for(i=0;ibegin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); - GetNodalFunctionsVelocity(vbasis, gauss); - - NormalBase(&bed_normal[0],&xyz_list_tria[0][0]); - bed_input->GetInputValue(&bed, gauss); - if(shelf_dampening){ //add dampening to avoid too high vertical velocities when not in hydrostatic equilibrium - bed_input->GetInputDerivativeValue(&dz[0],&xyz_list[0][0],gauss); - vx_input->GetInputValue(&vx, gauss); - vy_input->GetInputValue(&vy, gauss); - vz_input->GetInputValue(&vz, gauss); - dt=0.; - normal_vel=bed_normal[0]*vx+bed_normal[1]*vy+bed_normal[2]*vz; - damper=gravity*rho_water*pow(1+pow(dz[0],2)+pow(dz[1],2),0.5)*normal_vel*dt; - } - else damper=0.; - water_pressure=gravity*rho_water*bed; - - for(i=0;ivalues[3*i+0]+=(water_pressure+damper)*gauss->weight*Jdet2d*vbasis[i]*bed_normal[0]; - pe->values[3*i+1]+=(water_pressure+damper)*gauss->weight*Jdet2d*vbasis[i]*bed_normal[1]; - pe->values[3*i+2]+=(water_pressure+damper)*gauss->weight*Jdet2d*vbasis[i]*bed_normal[2]; - } - } - - /*Transform coordinate system*/ - ::TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list); - - /*Clean up and return*/ - delete gauss; - xDelete(cs_list); - xDelete(vbasis); - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorStressbalanceVert {{{*/ -ElementVector* Penta::CreatePVectorStressbalanceVert(void){ - - /*compute all load vectors for this element*/ - ElementVector* pe1=CreatePVectorStressbalanceVertVolume(); - ElementVector* pe2=CreatePVectorStressbalanceVertBase(); - ElementVector* pe =new ElementVector(pe1,pe2); - - /*clean-up and return*/ - delete pe1; - delete pe2; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorStressbalanceVertVolume {{{*/ -ElementVector* Penta::CreatePVectorStressbalanceVertVolume(void){ - - /*Constants*/ - const int numdof=NDOF1*NUMVERTICES; - - /*Intermediaries*/ - int approximation; - IssmDouble Jdet; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble dudx,dvdy,dwdz; - IssmDouble du[3],dv[3],dw[3]; - IssmDouble basis[6]; - GaussPenta *gauss=NULL; - - /*Initialize Element vector*/ - ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - inputs->GetInputValue(&approximation,ApproximationEnum); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vzFS_input=NULL; - if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){ - vzFS_input=inputs->GetInput(VzFSEnum); _assert_(vzFS_input); - } - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(2,2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctionsP1(basis, gauss); - - vx_input->GetInputDerivativeValue(&du[0],&xyz_list[0][0],gauss); - vy_input->GetInputDerivativeValue(&dv[0],&xyz_list[0][0],gauss); - if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){ - vzFS_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss); - dwdz=dw[2]; - } - else dwdz=0; - dudx=du[0]; - dvdy=dv[1]; - - for(int i=0;ivalues[i] += (dudx+dvdy+dwdz)*Jdet*gauss->weight*basis[i]; - } - - /*Clean up and return*/ - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorStressbalanceVertBase {{{*/ -ElementVector* Penta::CreatePVectorStressbalanceVertBase(void){ - - /*Constants*/ - const int numdof=NDOF1*NUMVERTICES; - - /*Intermediaries */ - int i,j; - int approximation; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_tria[NUMVERTICES2D][3]; - IssmDouble Jdet2d; - IssmDouble vx,vy,vz,dbdx,dbdy,basalmeltingvalue; - IssmDouble slope[3]; - IssmDouble basis[NUMVERTICES]; - GaussPenta* gauss=NULL; - - if (!IsOnBed()) return NULL; - - /*Initialize Element vector*/ - ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - for(i=0;iGetInputValue(&approximation,ApproximationEnum); - Input* bed_input=inputs->GetInput(BedEnum); _assert_(bed_input); - Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vzFS_input=NULL; - if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){ - vzFS_input=inputs->GetInput(VzFSEnum); _assert_(vzFS_input); - } - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(0,1,2,2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - basal_melting_input->GetInputValue(&basalmeltingvalue, gauss); - bed_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss); - vx_input->GetInputValue(&vx, gauss); - vy_input->GetInputValue(&vy, gauss); - if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){ - vzFS_input->GetInputValue(&vz, gauss); - } - else vz=0; - - dbdx=slope[0]; - dbdy=slope[1]; - - GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss); - GetNodalFunctionsP1(&basis[0], gauss); - - for(i=0;ivalues[i]+=-Jdet2d*gauss->weight*(vx*dbdx+vy*dbdy-vz-basalmeltingvalue)*basis[i]; - } - - /*Clean up and return*/ - delete gauss; - return pe; -} -/*}}}*/ /*FUNCTION Penta::CreateJacobianStressbalanceHoriz{{{*/ ElementMatrix* Penta::CreateJacobianStressbalanceHoriz(void){ @@ -9246,26 +7034,4 @@ } /*}}}*/ -/*FUNCTION Penta::CreatePVectorBalancethickness {{{*/ -ElementVector* Penta::CreatePVectorBalancethickness(void){ - - if (!IsOnBed()) return NULL; - - /*Depth Averaging Vx and Vy*/ - this->InputDepthAverageAtBase(VxEnum,VxAverageEnum); - this->InputDepthAverageAtBase(VyEnum,VyAverageEnum); - - /*Call Tria function*/ - Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. - ElementVector* pe=tria->CreatePVectorBalancethickness(); - delete tria->material; delete tria; - - /*Delete Vx and Vy averaged*/ - this->inputs->DeleteInput(VxAverageEnum); - this->inputs->DeleteInput(VyAverageEnum); - - /*Clean up and return*/ - return pe; -} -/*}}}*/ #endif @@ -9308,46 +7074,4 @@ /*clean up and return*/ return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorHydrologyDCInefficient {{{*/ -ElementVector* Penta::CreatePVectorHydrologyDCInefficient(void){ - - if (!IsOnBed()) return NULL; - - /*Call Tria function*/ - Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. - ElementVector* pe=tria->CreatePVectorHydrologyDCInefficient(); - delete tria->material; delete tria; - - /*Clean up and return*/ - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorHydrologyDCEfficient {{{*/ -ElementVector* Penta::CreatePVectorHydrologyDCEfficient(void){ - - if (!IsOnBed()) return NULL; - - /*Call Tria function*/ - Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. - ElementVector* pe=tria->CreatePVectorHydrologyDCEfficient(); - delete tria->material; delete tria; - - /*Clean up and return*/ - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorL@ProjectionEPL {{{*/ -ElementVector* Penta::CreatePVectorL2ProjectionEPL(void){ - - if (!IsOnBed()) return NULL; - - /*Call Tria function*/ - Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. - ElementVector* pe=tria->CreatePVectorL2ProjectionEPL(); - delete tria->material; delete tria; - - /*Clean up and return*/ - return pe; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16981) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16982) @@ -199,9 +199,4 @@ ElementMatrix* CreateKMatrixFreeSurfaceTop(void); ElementMatrix* CreateKMatrixFreeSurfaceBase(void); - ElementVector* CreatePVector(void); - ElementVector* CreatePVectorMasstransport(void); - ElementVector* CreatePVectorFreeSurfaceTop(void); - ElementVector* CreatePVectorFreeSurfaceBase(void); - ElementVector* CreatePVectorL2ProjectionBase(void); IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure); IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure); @@ -301,38 +296,10 @@ ElementMatrix* CreateJacobianStressbalanceHO(void); ElementMatrix* CreateJacobianStressbalanceFS(void); - ElementVector* CreatePVectorCouplingSSAFS(void); - ElementVector* CreatePVectorCouplingSSAFSViscous(void); - ElementVector* CreatePVectorCouplingSSAFSFriction(void); - ElementVector* CreatePVectorCouplingHOFS(void); - ElementVector* CreatePVectorCouplingHOFSViscous(void); - ElementVector* CreatePVectorCouplingHOFSFriction(void); - ElementVector* CreatePVectorStressbalanceHoriz(void); - ElementVector* CreatePVectorStressbalanceSIA(void); - ElementVector* CreatePVectorStressbalanceSSA(void); - ElementVector* CreatePVectorStressbalanceSSAHO(void); - ElementVector* CreatePVectorStressbalanceSSAFS(void); - ElementVector* CreatePVectorStressbalanceL1L2(void); - ElementVector* CreatePVectorStressbalanceHO(void); - ElementVector* CreatePVectorStressbalanceHODrivingStress(void); - ElementVector* CreatePVectorStressbalanceHOFront(void); - ElementVector* CreatePVectorStressbalanceHOFS(void); - ElementVector* CreatePVectorStressbalanceFS(void); - ElementVector* CreatePVectorStressbalanceFSFront(void); - ElementVector* CreatePVectorStressbalanceFSViscous(void); - void PVectorGLSstabilization(ElementVector* pe); - ElementVector* CreatePVectorStressbalanceFSShelf(void); - ElementVector* CreatePVectorStressbalanceVert(void); - ElementVector* CreatePVectorStressbalanceVertVolume(void); - ElementVector* CreatePVectorStressbalanceVertBase(void); #endif #ifdef _HAVE_CONTROL_ - ElementVector* CreatePVectorAdjointHoriz(void); ElementMatrix* CreateKMatrixAdjointSSA2d(void); ElementMatrix* CreateKMatrixAdjointHO(void); ElementMatrix* CreateKMatrixAdjointFS(void); - ElementVector* CreatePVectorAdjointSSA(void); - ElementVector* CreatePVectorAdjointHO(void); - ElementVector* CreatePVectorAdjointFS(void); #endif @@ -341,7 +308,4 @@ ElementMatrix* CreateKMatrixHydrologyDCEfficient(void); ElementMatrix* CreateEPLDomainMassMatrix(void); - ElementVector* CreatePVectorHydrologyDCInefficient(void); - ElementVector* CreatePVectorHydrologyDCEfficient(void); - ElementVector* CreatePVectorL2ProjectionEPL(void); void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode); void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index){_error_("not implemented yet");}; @@ -362,13 +326,4 @@ ElementMatrix* CreateKMatrixThermalVolume(void); ElementMatrix* CreateKMatrixThermalShelf(void); - ElementVector* CreatePVectorEnthalpy(void); - ElementVector* CreatePVectorEnthalpyVolume(void); - ElementVector* CreatePVectorEnthalpyShelf(void); - ElementVector* CreatePVectorEnthalpySheet(void); - ElementVector* CreatePVectorMelting(void); - ElementVector* CreatePVectorThermal(void); - ElementVector* CreatePVectorThermalVolume(void); - ElementVector* CreatePVectorThermalShelf(void); - ElementVector* CreatePVectorThermalSheet(void); void UpdateBasalConstraintsEnthalpy(void); void ComputeBasalMeltingrate(void); @@ -378,5 +333,4 @@ #ifdef _HAVE_BALANCED_ ElementMatrix* CreateKMatrixBalancethickness(void); - ElementVector* CreatePVectorBalancethickness(void); #endif /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp (revision 16981) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp (revision 16982) @@ -333,151 +333,4 @@ } /*}}}*/ -/*FUNCTION Seg::CreatePVectorL2Projection {{{*/ -ElementVector* Seg::CreatePVectorL2Projection(void){ - - /*Intermediaries */ - int i,input_enum; - IssmDouble Jdet; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble value; - Input* input = NULL; - Input* input2 = NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element vector*/ - ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&input_enum,InputToL2ProjectEnum); - switch(input_enum){ - case SurfaceSlopeXEnum: input2 = inputs->GetInput(SurfaceEnum); _assert_(input2); break; - case BedSlopeXEnum: input2 = inputs->GetInput(BedEnum); _assert_(input2); break; - default: input = inputs->GetInput(input_enum); - } - - /* Start looping on the number of gaussian points: */ - GaussSeg* gauss=new GaussSeg(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - - if(input2){ - input2->GetInputDerivativeValue(&value,&xyz_list[0][0],gauss); - } - else{ - input->GetInputValue(&value,gauss); - } - - for(i=0;ivalues[i]+=Jdet*gauss->weight*value*basis[i]; - } - - /*Clean up and return*/ - xDelete(basis); - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Seg::CreatePVectorFreeSurfaceTop {{{*/ -ElementVector* Seg::CreatePVectorFreeSurfaceTop(void){ - - /*Intermediaries */ - IssmDouble Jdet,dt; - IssmDouble ms,surface,vy; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element vector and other vectors*/ - ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input); - Input* ms_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(ms_input); - Input* surface_input= inputs->GetInput(SurfaceEnum); _assert_(surface_input); - - /*Initialize mb_correction to 0, do not forget!:*/ - /* Start looping on the number of gaussian points: */ - GaussSeg* gauss=new GaussSeg(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - - vy_input->GetInputValue(&vy,gauss); - ms_input->GetInputValue(&ms,gauss); - surface_input->GetInputValue(&surface,gauss); - - for(int i=0;ivalues[i]+=Jdet*gauss->weight*(surface + dt*ms + dt*vy)*basis[i]; - } - - /*Clean up and return*/ - xDelete(basis); - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Seg::CreatePVectorFreeSurfaceBase {{{*/ -ElementVector* Seg::CreatePVectorFreeSurfaceBase(void){ - - /*Intermediaries */ - IssmDouble Jdet,dt; - IssmDouble mb,mb_correction,bed,vy; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element vector and other vectors*/ - ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input); - Input* mb_input = inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(mb_input); - Input* mb_correction_input = inputs->GetInput(BasalforcingsMeltingRateCorrectionEnum); - Input* bed_input = inputs->GetInput(BedEnum); _assert_(bed_input); - - /*Initialize mb_correction to 0, do not forget!:*/ - /* Start looping on the number of gaussian points: */ - GaussSeg* gauss=new GaussSeg(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - - vy_input->GetInputValue(&vy,gauss); - mb_input->GetInputValue(&mb,gauss); - bed_input->GetInputValue(&bed,gauss); - if(mb_correction_input) - mb_correction_input->GetInputValue(&mb_correction,gauss); - else - mb_correction=0.; - - for(int i=0;ivalues[i]+=Jdet*gauss->weight*(bed+dt*(mb-mb_correction) + dt*vy)*basis[i]; - } - - /*Clean up and return*/ - xDelete(basis); - delete gauss; - return pe; -} -/*}}}*/ /*FUNCTION Seg::GetNumberOfNodes;{{{*/ int Seg::GetNumberOfNodes(void){ Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16981) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16982) @@ -71,6 +71,4 @@ ElementMatrix* CreateKMatrix(void){_error_("not implemented yet");}; void CreateDVector(Vector* df){_error_("not implemented yet");}; - ElementVector* CreatePVector(void){_error_("not implemented yet");}; - ElementVector* CreatePVectorL2Projection(void); void CreateJacobianMatrix(Matrix* Jff){_error_("not implemented yet");}; void Delta18oParameterization(void){_error_("not implemented yet");}; @@ -262,6 +260,4 @@ ElementMatrix* CreateKMatrixFreeSurfaceTop(void); ElementMatrix* CreateKMatrixFreeSurfaceBase(void); - ElementVector* CreatePVectorFreeSurfaceTop(void); - ElementVector* CreatePVectorFreeSurfaceBase(void); IssmDouble GetSize(void); }; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16981) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16982) @@ -374,201 +374,4 @@ } /*}}}*/ -/*FUNCTION Tria::CreatePVector(void){{{*/ -ElementVector* Tria::CreatePVector(void){ - - /*retrive parameters: */ - int meshtype; - int analysis_type; - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*asserts: {{{*/ - /*if debugging mode, check that all pointers exist*/ - _assert_(this->nodes && this->material && this->matpar && this->parameters && this->inputs); - /*}}}*/ - - /*Skip if water element*/ - if(NoIceInElement()) return NULL; - - /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */ - switch(analysis_type){ - #ifdef _HAVE_STRESSBALANCE_ - case StressbalanceAnalysisEnum: - int approximation; - inputs->GetInputValue(&approximation,ApproximationEnum); - switch(approximation){ - case SSAApproximationEnum: - return CreatePVectorStressbalanceSSA(); - case FSApproximationEnum: - return CreatePVectorStressbalanceFS(); - case SIAApproximationEnum: - return NULL; - case NoneApproximationEnum: - return NULL; - default: - _error_("Approximation " << EnumToStringx(approximation) << " not supported yet"); - } - break; - case StressbalanceSIAAnalysisEnum: - return CreatePVectorStressbalanceSIA(); - break; - #endif - #ifdef _HAVE_DAMAGE_ - case DamageEvolutionAnalysisEnum: - return CreatePVectorDamageEvolution(); - break; - #endif - case L2ProjectionBaseAnalysisEnum: - parameters->FindParam(&meshtype,MeshTypeEnum); - if(meshtype==Mesh2DverticalEnum){ - return CreatePVectorL2ProjectionBase(); - } - else{ - return CreatePVectorL2Projection(); - } - break; - #ifdef _HAVE_MASSTRANSPORT_ - case MasstransportAnalysisEnum: - return CreatePVectorMasstransport(); - break; - case FreeSurfaceTopAnalysisEnum: - parameters->FindParam(&meshtype,MeshTypeEnum); - switch(meshtype){ - case Mesh2DverticalEnum: - return CreatePVectorFreeSurfaceTop1D(); - case Mesh3DEnum: - return CreatePVectorFreeSurfaceTop(); - default: - _error_("Mesh not supported yet"); - } - break; - case FreeSurfaceBaseAnalysisEnum: - parameters->FindParam(&meshtype,MeshTypeEnum); - switch(meshtype){ - case Mesh2DverticalEnum: - return CreatePVectorFreeSurfaceBase1D(); - case Mesh3DEnum: - return CreatePVectorFreeSurfaceBase(); - default: - _error_("Mesh not supported yet"); - } - break; - case ExtrudeFromBaseAnalysisEnum: case ExtrudeFromTopAnalysisEnum: - return NULL; - break; - #endif - #ifdef _HAVE_HYDROLOGY_ - case HydrologyShreveAnalysisEnum: - return CreatePVectorHydrologyShreve(); - break; - case HydrologyDCInefficientAnalysisEnum: - return CreatePVectorHydrologyDCInefficient(); - break; - case HydrologyDCEfficientAnalysisEnum: - return CreatePVectorHydrologyDCEfficient(); - break; - case L2ProjectionEPLAnalysisEnum: - return CreatePVectorL2ProjectionEPL(); - break; - #endif - #ifdef _HAVE_BALANCED_ - case BalancethicknessAnalysisEnum: - return CreatePVectorBalancethickness(); - break; - case BalancevelocityAnalysisEnum: - return CreatePVectorBalancevelocity(); - break; - case SmoothedSurfaceSlopeXAnalysisEnum: - return CreatePVectorSmoothedSlopeX(); - break; - case SmoothedSurfaceSlopeYAnalysisEnum: - return CreatePVectorSmoothedSlopeY(); - break; - #endif - #ifdef _HAVE_CONTROL_ - case AdjointBalancethicknessAnalysisEnum: - return CreatePVectorAdjointBalancethickness(); - break; - case AdjointHorizAnalysisEnum: - return CreatePVectorAdjointHoriz(); - break; - #endif - default: - _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet"); - } - - /*make compiler happy*/ - return NULL; -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorL2ProjectionBase{{{*/ -ElementVector* Tria::CreatePVectorL2ProjectionBase(void){ - - if(!HasEdgeOnBed()) return NULL; - - int index1,index2; - this->EdgeOnBedIndices(&index1,&index2); - - Seg* seg=(Seg*)SpawnSeg(index1,index2); - ElementVector* pe=seg->CreatePVectorL2Projection(); - delete seg->material; delete seg; - - /*clean up and return*/ - return pe; -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorL2Projection {{{*/ -ElementVector* Tria::CreatePVectorL2Projection(void){ - - /*Intermediaries */ - int i,input_enum; - IssmDouble Jdet,value; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble slopes[2]; - Input* input = NULL; - Input* input2 = NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element vector*/ - ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&input_enum,InputToL2ProjectEnum); - switch(input_enum){ - case SurfaceSlopeXEnum: input2 = inputs->GetInput(SurfaceEnum); _assert_(input2); break; - case SurfaceSlopeYEnum: input2 = inputs->GetInput(SurfaceEnum); _assert_(input2); break; - case BedSlopeXEnum: input2 = inputs->GetInput(BedEnum); _assert_(input2); break; - case BedSlopeYEnum: input2 = inputs->GetInput(BedEnum); _assert_(input2); break; - default: input = inputs->GetInput(input_enum); - } - - /* Start looping on the number of gaussian points: */ - GaussTria* gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - - if(input2) input2->GetInputDerivativeValue(&slopes[0],&xyz_list[0][0],gauss); - switch(input_enum){ - case SurfaceSlopeXEnum: case BedSlopeXEnum: value = slopes[0]; break; - case SurfaceSlopeYEnum: case BedSlopeYEnum: value = slopes[1]; break; - default: input->GetInputValue(&value,gauss); - } - - for(i=0;ivalues[i]+=Jdet*gauss->weight*value*basis[i]; - } - - /*Clean up and return*/ - xDelete(basis); - delete gauss; - return pe; -} -/*}}}*/ /*FUNCTION Tria::CreateJacobianMatrix{{{*/ void Tria::CreateJacobianMatrix(Matrix* Jff){ @@ -3817,420 +3620,4 @@ } /*}}}*/ -/*FUNCTION Tria::CreatePVectorStressbalanceFS {{{*/ -ElementVector* Tria::CreatePVectorStressbalanceFS(void){ - - ElementVector* pe1; - ElementVector* pe2; - ElementVector* pe3; - ElementVector* pe; - - /*compute all stiffness matrices for this element*/ - pe1=CreatePVectorStressbalanceFSViscous(); - pe2=CreatePVectorStressbalanceFSShelf(); - pe3=CreatePVectorStressbalanceFSFront(); - pe =new ElementVector(pe1,pe2,pe3); - - /*clean-up and return*/ - delete pe1; - delete pe2; - delete pe3; - return pe; -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorStressbalanceFSFront{{{*/ -ElementVector* Tria::CreatePVectorStressbalanceFSFront(void){ - - /*If no front, return NULL*/ - if(!IsZeroLevelset(MaskIceLevelsetEnum)) return NULL; - - /*Intermediaries*/ - int i; - IssmDouble rho_ice,rho_water,gravity,y_g; - IssmDouble Jdet,water_pressure,air_pressure,pressure; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble* xyz_list_front = NULL; - IssmDouble area_coordinates[2][3]; - IssmDouble normal[2]; - GaussTria* gauss = NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = this->NumberofNodesVelocity(); - int pnumnodes = this->NumberofNodesPressure(); - - /*Prepare coordinate system list*/ - int* cs_list = xNew(vnumnodes+pnumnodes); - for(i=0;iparameters,FSvelocityEnum); - IssmDouble* vbasis = xNew(vnumnodes); - - /*Retrieve all inputs and parameters*/ - rho_water = matpar->GetRhoWater(); - rho_ice = matpar->GetRhoIce(); - gravity = matpar->GetG(); - ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum); - GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],2); - NormalSection(&normal[0],xyz_list_front); - - /*Start looping on Gaussian points*/ - IssmDouble ymax=max(xyz_list_front[0*3+1],xyz_list_front[1*3+1]); - IssmDouble ymin=min(xyz_list_front[0*3+1],xyz_list_front[1*3+1]); - if(ymax>0. && ymin<0.) gauss=new GaussTria(area_coordinates,30); //refined in vertical because of the sea level discontinuity - else gauss=new GaussTria(area_coordinates,3); - - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - y_g=GetYcoord(gauss); - GetNodalFunctionsVelocity(vbasis,gauss); - GetSegmentJacobianDeterminant(&Jdet,xyz_list_front,gauss); - - water_pressure = rho_water*gravity*min(0.,y_g); //0 if the gaussian point is above water level - air_pressure = 0.; - pressure = water_pressure + air_pressure; - - for(i=0;ivalues[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*vbasis[i]; - pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*vbasis[i]; - } - } - - /*Transform coordinate system*/ - ::TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list); - - /*Clean up and return*/ - xDelete(cs_list); - xDelete(xyz_list_front); - xDelete(vbasis); - delete gauss; - return pe; - -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorStressbalanceFSViscous {{{*/ -ElementVector* Tria::CreatePVectorStressbalanceFSViscous(void){ - - /*Intermediaries*/ - int i; - int approximation; - IssmDouble Jdet,gravity,rho_ice; - IssmDouble forcex,forcey; - IssmDouble xyz_list[NUMVERTICES][3]; - GaussTria *gauss=NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = this->NumberofNodesVelocity(); - int pnumnodes = this->NumberofNodesPressure(); - - /*Prepare coordinate system list*/ - int* cs_list = xNew(vnumnodes+pnumnodes); - for(i=0;iparameters,FSvelocityEnum); - IssmDouble* vbasis = xNew(vnumnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* loadingforcex_input=inputs->GetInput(LoadingforceXEnum); _assert_(loadingforcex_input); - Input* loadingforcey_input=inputs->GetInput(LoadingforceYEnum); _assert_(loadingforcey_input); - rho_ice = matpar->GetRhoIce(); - gravity = matpar->GetG(); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(5); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctionsVelocity(vbasis, gauss); - - loadingforcex_input->GetInputValue(&forcex,gauss); - loadingforcey_input->GetInputValue(&forcey,gauss); - - for(i=0;ivalues[i*NDOF2+1] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i]; - pe->values[i*NDOF2+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i]; - pe->values[i*NDOF2+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i]; - } - } - - /*Transform coordinate system*/ - ::TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list); - - /*Clean up and return*/ - delete gauss; - xDelete(cs_list); - xDelete(vbasis); - return pe; -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorStressbalanceFSShelf{{{*/ -ElementVector* Tria::CreatePVectorStressbalanceFSShelf(void){ - - /*Intermediaries*/ - int i,j; - int indices[2]; - IssmDouble gravity,rho_water,bed,water_pressure; - IssmDouble normal_vel,vx,vy,dt; - IssmDouble xyz_list_seg[NUMVERTICES1D][3]; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble normal[2]; - IssmDouble Jdet; - - /*Initialize Element vector and return if necessary*/ - if(!HasEdgeOnBed() || !IsFloating()) return NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = this->NumberofNodesVelocity(); - int pnumnodes = this->NumberofNodesPressure(); - - /*Prepare coordinate system list*/ - int* cs_list = xNew(vnumnodes+pnumnodes); - for(i=0;iparameters,FSvelocityEnum); - IssmDouble* vbasis = xNew(vnumnodes); - - /*Retrieve all inputs and parameters*/ - rho_water=matpar->GetRhoWater(); - gravity=matpar->GetG(); - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* bed_input=inputs->GetInput(BedEnum); _assert_(bed_input); - - /*Get vertex indices that lie on bed*/ - this->EdgeOnBedIndices(&indices[0],&indices[1]); - - for(i=0;ibegin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetSegmentJacobianDeterminant(&Jdet,&xyz_list_seg[0][0],gauss); - GetNodalFunctionsVelocity(vbasis, gauss); - - NormalSection(&normal[0],&xyz_list_seg[0][0]); - _assert_(normal[1]<0.); - bed_input->GetInputValue(&bed, gauss); - water_pressure=gravity*rho_water*bed; - - for(i=0;ivalues[2*i+0]+=(water_pressure)*gauss->weight*Jdet*vbasis[i]*normal[0]; - pe->values[2*i+1]+=(water_pressure)*gauss->weight*Jdet*vbasis[i]*normal[1]; - } - } - - /*Transform coordinate system*/ - ::TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list); - - /*Clean up and return*/ - delete gauss; - xDelete(cs_list); - xDelete(vbasis); - return pe; -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorStressbalanceSSA {{{*/ -ElementVector* Tria::CreatePVectorStressbalanceSSA(){ - - /*compute all load vectors for this element*/ - ElementVector* pe1=CreatePVectorStressbalanceSSADrivingStress(); - ElementVector* pe2=CreatePVectorStressbalanceSSAFront(); - ElementVector* pe =new ElementVector(pe1,pe2); - - /*clean-up and return*/ - delete pe1; - delete pe2; - return pe; -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorStressbalanceSSADrivingStress {{{*/ -ElementVector* Tria::CreatePVectorStressbalanceSSADrivingStress(){ - - /*Intermediaries */ - int i; - IssmDouble driving_stress_baseline,thickness; - IssmDouble Jdet; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble slope[2]; - IssmDouble icefrontlevel[3]; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element vector and vectors*/ - ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters,SSAApproximationEnum); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); - Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input); - Input* drag_input=inputs->GetInput(FrictionCoefficientEnum);_assert_(drag_input); - GetInputListOnVertices(&icefrontlevel[0],BedEnum); - - /* Start looping on the number of gaussian points: */ - GaussTria* gauss = new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis, gauss); - - thickness_input->GetInputValue(&thickness,gauss); - surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss); - driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG()*thickness; - - /*Build load vector: */ - for(i=0;ivalues[i*NDOF2+0]+=-driving_stress_baseline*slope[0]*Jdet*gauss->weight*basis[i]; - pe->values[i*NDOF2+1]+=-driving_stress_baseline*slope[1]*Jdet*gauss->weight*basis[i]; - } - } - - /*Transform coordinate system*/ - ::TransformLoadVectorCoord(pe,nodes,numnodes,XYEnum); - - /*Clean up and return*/ - xDelete(basis); - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorStressbalanceSSAFront {{{*/ -ElementVector* Tria::CreatePVectorStressbalanceSSAFront(){ - - /*If no front, return NULL*/ - if(!IsZeroLevelset(MaskIceLevelsetEnum)) return NULL; - - /*Intermediaries*/ - IssmDouble rho_ice,rho_water,gravity; - IssmDouble Jdet,thickness,bed,water_pressure,ice_pressure,air_pressure; - IssmDouble surface_under_water,base_under_water,pressure; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble *xyz_list_front = NULL; - IssmDouble area_coordinates[2][3]; - IssmDouble normal[2]; - - /*Fetch number of nodes for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element vector and other vectors*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters,SSAApproximationEnum); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); - Input* bed_input =inputs->GetInput(BedEnum); _assert_(bed_input); - rho_water = matpar->GetRhoWater(); - rho_ice = matpar->GetRhoIce(); - gravity = matpar->GetG(); - ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum); - GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],2); - NormalSection(&normal[0],xyz_list_front); - - /*Start looping on Gaussian points*/ - GaussTria* gauss=new GaussTria(area_coordinates,3); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - thickness_input->GetInputValue(&thickness,gauss); - bed_input->GetInputValue(&bed,gauss); - GetSegmentJacobianDeterminant(&Jdet,xyz_list_front,gauss); - GetNodalFunctions(basis,gauss); - - surface_under_water = min(0.,thickness+bed); // 0 if the top of the glacier is above water level - base_under_water = min(0.,bed); // 0 if the bottom of the glacier is above water level - water_pressure = 1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2)); - ice_pressure = 1.0/2.0*gravity*rho_ice*pow(thickness,2); - air_pressure = 0; - pressure = ice_pressure + water_pressure + air_pressure; - - for (int i=0;ivalues[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i]; - pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i]; - } - } - - /*Transform coordinate system*/ - ::TransformLoadVectorCoord(pe,nodes,numnodes,XYEnum); - - /*Clean up and return*/ - xDelete(basis); - xDelete(xyz_list_front); - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorStressbalanceSIA{{{*/ -ElementVector* Tria::CreatePVectorStressbalanceSIA(void){ - - /*Intermediaries */ - IssmDouble ub,vb; - IssmDouble rho_ice,gravity,n,B,drag; - IssmDouble slope2,thickness,connectivity; - IssmDouble slope[2]; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); _assert_(numnodes==3); - - /*Initialize Element vector*/ - ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters); - - /*Retrieve all inputs and parameters*/ - rho_ice=matpar->GetRhoIce(); - gravity=matpar->GetG(); - n=material->GetN(); - B=material->GetBbar(); - Input* slopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(slopex_input); - Input* slopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(slopey_input); - Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); - Input* drag_input=inputs->GetInput(FrictionCoefficientEnum); _assert_(drag_input); - - /*Spawn 3 sing elements: */ - GaussTria* gauss=new GaussTria(); - for(int i=0;iGaussVertex(i); - - connectivity=(IssmDouble)vertices[i]->Connectivity(); - - thickness_input->GetInputValue(&thickness,gauss); - drag_input->GetInputValue(&drag,gauss); - slopex_input->GetInputValue(&slope[0],gauss); - slopey_input->GetInputValue(&slope[1],gauss); - slope2=slope[0]*slope[0]+slope[1]*slope[1]; - - /*Payne 1995 (m = 1, B = 5e-3/yts = 1.58e-10 m/s/Pa*/ - ub=-1.58*1.e-10*rho_ice*gravity*thickness*slope[0]; - vb=-1.58*1.e-10*rho_ice*gravity*thickness*slope[1]; - ///*Ritz et al. 1996*/ - //ub=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[0]/sqrt(slope2); - //vb=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[1]/sqrt(slope2); - ///*Rutt et al. 2009*/ - //ub=-drag*rho_ice*gravity*thickness*slope[0]; - //vb=-drag*rho_ice*gravity*thickness*slope[1]; - - pe->values[2*i+0]=(ub-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[0])/connectivity; - pe->values[2*i+1]=(vb-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[1])/connectivity; - } - - /*Clean up and return*/ - delete gauss; - return pe; -} -/*}}}*/ /*FUNCTION Tria::CreateJacobianStressbalanceSSA{{{*/ ElementMatrix* Tria::CreateJacobianStressbalanceSSA(void){ @@ -5607,277 +4994,4 @@ } /*}}}*/ -/*FUNCTION Tria::CreatePVectorAdjointBalancethickness{{{*/ -ElementVector* Tria::CreatePVectorAdjointBalancethickness(void){ - - /*Intermediaries */ - int i,resp; - IssmDouble Jdet; - IssmDouble thickness,thicknessobs,weight; - int num_responses; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble dH[2]; - IssmDouble vx,vy,vel; - int *responses = NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element vector and vectors*/ - ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); - IssmDouble* basis = xNew(numnodes); - IssmDouble* dbasis = xNew(2*numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum); - this->parameters->FindParam(&responses,NULL,InversionCostFunctionsEnum); - Input* thickness_input = inputs->GetInput(ThicknessEnum); _assert_(thickness_input); - Input* thicknessobs_input = inputs->GetInput(InversionThicknessObsEnum); _assert_(thicknessobs_input); - Input* weights_input = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); - Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input); - - /* Start looping on the number of gaussian points: */ - GaussTria* gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis, gauss); - GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss); - - thickness_input->GetInputValue(&thickness, gauss); - thickness_input->GetInputDerivativeValue(&dH[0],&xyz_list[0][0],gauss); - thicknessobs_input->GetInputValue(&thicknessobs, gauss); - - /*Loop over all requested responses*/ - for(resp=0;respGetInputValue(&weight,gauss,responses[resp]); - - switch(responses[resp]){ - case ThicknessAbsMisfitEnum: - for(i=0;ivalues[i]+=(thicknessobs-thickness)*weight*Jdet*gauss->weight*basis[i]; - break; - case ThicknessAbsGradientEnum: - for(i=0;ivalues[i]+= - weight*dH[0]*dbasis[0*numnodes+i]*Jdet*gauss->weight; - for(i=0;ivalues[i]+= - weight*dH[1]*dbasis[1*numnodes+i]*Jdet*gauss->weight; - break; - case ThicknessAlongGradientEnum: - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vel = sqrt(vx*vx+vy*vy); - vx = vx/(vel+1.e-9); - vy = vy/(vel+1.e-9); - for(i=0;ivalues[i]+= - weight*(dH[0]*vx+dH[1]*vy)*(dbasis[0*numnodes+i]*vx+dbasis[1*numnodes+i]*vy)*Jdet*gauss->weight; - break; - case ThicknessAcrossGradientEnum: - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vel = sqrt(vx*vx+vy*vy); - vx = vx/(vel+1.e-9); - vy = vy/(vel+1.e-9); - for(i=0;ivalues[i]+= - weight*(dH[0]*(-vy)+dH[1]*vx)*(dbasis[0*numnodes+i]*(-vy)+dbasis[1*numnodes+i]*vx)*Jdet*gauss->weight; - break; - default: - _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); - } - } - } - - /*Clean up and return*/ - xDelete(basis); - xDelete(dbasis); - xDelete(responses); - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorAdjointHoriz{{{*/ -ElementVector* Tria::CreatePVectorAdjointHoriz(void){ - - /*Intermediaries */ - int i,resp; - int *responses=NULL; - int num_responses; - IssmDouble Jdet; - IssmDouble obs_velocity_mag,velocity_mag; - IssmDouble dux,duy; - IssmDouble epsvel=2.220446049250313e-16; - IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/ - IssmDouble scalex=0.,scaley=0.,scale=0.,S=0.; - IssmDouble vx,vy,vxobs,vyobs,weight; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element vector*/ - ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum); - this->parameters->FindParam(&responses,NULL,InversionCostFunctionsEnum); - Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); - Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input =inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vxobs_input =inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input); - Input* vyobs_input =inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input); - - /*Get Surface if required by one response*/ - for(resp=0;respGetInputValue(&S,SurfaceAreaEnum); break; - } - } - - /* Start looping on the number of gaussian points: */ - GaussTria* gauss=new GaussTria(4); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - - /*Get all parameters at gaussian point*/ - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vxobs_input->GetInputValue(&vxobs,gauss); - vyobs_input->GetInputValue(&vyobs,gauss); - GetNodalFunctions(basis, gauss); - - /*Loop over all requested responses*/ - for(resp=0;respGetInputValue(&weight,gauss,responses[resp]); - - switch(responses[resp]){ - case SurfaceAbsVelMisfitEnum: - /* - * 1 [ 2 2 ] - * J = --- | (u - u ) + (v - v ) | - * 2 [ obs obs ] - * - * dJ - * DU = - -- = (u - u ) - * du obs - */ - for(i=0;ivalues[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; - pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; - } - break; - case SurfaceRelVelMisfitEnum: - /* - * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] - * J = --- | ------------- (u - u ) + ------------- (v - v ) | - * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] - * obs obs - * - * dJ \bar{v}^2 - * DU = - -- = ------------- (u - u ) - * du (u + eps)^2 obs - * obs - */ - for(i=0;ivalues[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; - pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; - } - break; - case SurfaceLogVelMisfitEnum: - /* - * [ vel + eps ] 2 - * J = 4 \bar{v}^2 | log ( ----------- ) | - * [ vel + eps ] - * obs - * - * dJ 2 * log(...) - * DU = - -- = - 4 \bar{v}^2 ------------- u - * du vel^2 + eps - * - */ - for(i=0;ivalues[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; - pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; - } - break; - case SurfaceAverageVelMisfitEnum: - /* - * 1 2 2 - * J = --- sqrt( (u - u ) + (v - v ) ) - * S obs obs - * - * dJ 1 1 - * DU = - -- = - --- ----------- * 2 (u - u ) - * du S 2 sqrt(...) obs - */ - for(i=0;ivalues[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; - pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; - } - break; - case SurfaceLogVxVyMisfitEnum: - /* - * 1 [ |u| + eps 2 |v| + eps 2 ] - * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | - * 2 [ |u |+ eps |v |+ eps ] - * obs obs - * dJ 1 u 1 - * DU = - -- = - \bar{v}^2 log(u...) --------- ---- ~ - \bar{v}^2 log(u...) ------ - * du |u| + eps |u| u + eps - */ - for(i=0;ivalues[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; - pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; - } - break; - case DragCoefficientAbsGradientEnum: - /*Nothing in P vector*/ - break; - case ThicknessAbsGradientEnum: - /*Nothing in P vector*/ - break; - case ThicknessAlongGradientEnum: - /*Nothing in P vector*/ - break; - case ThicknessAcrossGradientEnum: - /*Nothing in P vector*/ - break; - case RheologyBbarAbsGradientEnum: - /*Nothing in P vector*/ - break; - default: - _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); - } - } - } - - /*Clean up and return*/ - xDelete(basis); - xDelete(responses); - delete gauss; - return pe; -} -/*}}}*/ /*FUNCTION Tria::DragCoefficientAbsGradient{{{*/ IssmDouble Tria::DragCoefficientAbsGradient(void){ @@ -6449,244 +5563,4 @@ xDelete(basis); return Ke; -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorHydrologyShreve {{{*/ -ElementVector* Tria::CreatePVectorHydrologyShreve(void){ - - /*Intermediaries */ - int i; - IssmDouble Jdet,dt; - IssmDouble basal_melting_g; - IssmDouble old_watercolumn_g; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Skip if water or ice shelf element*/ - if(NoIceInElement() || IsFloating()) return NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element vector*/ - ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters); - IssmDouble* basis = xNewZeroInit(numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input); - Input* old_watercolumn_input=inputs->GetInput(WaterColumnOldEnum); _assert_(old_watercolumn_input); - /*Initialize basal_melting_correction_g to 0, do not forget!:*/ - /* Start looping on the number of gaussian points: */ - GaussTria* gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis, gauss); - - basal_melting_input->GetInputValue(&basal_melting_g,gauss); - old_watercolumn_input->GetInputValue(&old_watercolumn_g,gauss); - - if(reCast(dt)){ - for(i=0;ivalues[i]+=Jdet*gauss->weight*(old_watercolumn_g+dt*basal_melting_g)*basis[i]; - } - else{ - for(i=0;ivalues[i]+=Jdet*gauss->weight*basal_melting_g*basis[i]; - } - } - - /*Clean up and return*/ - xDelete(basis); - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorHydrologyDCInefficient {{{*/ -ElementVector* Tria::CreatePVectorHydrologyDCInefficient(void){ - - /*Intermediaries */ - IssmDouble Jdet; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble dt,scalar,water_head; - IssmDouble water_load,transfer; - IssmDouble sediment_storing; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element vector and other vectors*/ - ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - sediment_storing = matpar->GetSedimentStoring(); - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - Input* water_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(water_input); - Input* transfer_input=inputs->GetInput(WaterTransferEnum); _assert_(transfer_input); - Input* old_wh_input=NULL; - - if(reCast(dt)){ - old_wh_input=inputs->GetInput(SedimentHeadOldEnum); _assert_(old_wh_input); - } - - /* Start looping on the number of gaussian points: */ - GaussTria* gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis, gauss); - - /*Loading term*/ - water_input->GetInputValue(&water_load,gauss); - transfer_input->GetInputValue(&transfer,gauss); - scalar = Jdet*gauss->weight*(water_load+transfer); - if(reCast(dt)) scalar = scalar*dt; - for(int i=0;ivalues[i]+=scalar*basis[i]; - - /*Transient term*/ - if(reCast(dt)){ - old_wh_input->GetInputValue(&water_head,gauss); - scalar = Jdet*gauss->weight*water_head*sediment_storing; - for(int i=0;ivalues[i]+=scalar*basis[i]; - } - } - - /*Clean up and return*/ - xDelete(basis); - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorHydrologyDCEfficient {{{*/ -ElementVector* Tria::CreatePVectorHydrologyDCEfficient(void){ - - /*Intermediaries */ - IssmDouble connectivity; - IssmDouble Jdet; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble dt,scalar,water_head; - IssmDouble transfer,residual; - IssmDouble epl_thickness; - IssmDouble epl_specificstoring; - GaussTria* gauss = NULL; - - /*Check that all nodes are active, else return empty matrix*/ - if(!this->AllActive()){ - return NULL; - } - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element vector and other vectors*/ - ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - epl_specificstoring = matpar->GetEplSpecificStoring(); - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - Input* residual_input=inputs->GetInput(SedimentHeadResidualEnum); _assert_(residual_input); - Input* transfer_input=inputs->GetInput(WaterTransferEnum); _assert_(transfer_input); - Input* thickness_input=inputs->GetInput(HydrologydcEplThicknessEnum); _assert_(thickness_input); - Input* old_wh_input=NULL; - - if(reCast(dt)){ - old_wh_input=inputs->GetInput(EplHeadOldEnum); _assert_(old_wh_input); - } - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - - /*Loading term*/ - transfer_input->GetInputValue(&transfer,gauss); - scalar = Jdet*gauss->weight*(-transfer); - if(reCast(dt)) scalar = scalar*dt; - for(int i=0;ivalues[i]+=scalar*basis[i]; - - /*Transient term*/ - if(reCast(dt)){ - thickness_input->GetInputValue(&epl_thickness,gauss); - old_wh_input->GetInputValue(&water_head,gauss); - scalar = Jdet*gauss->weight*water_head*epl_specificstoring*epl_thickness; - for(int i=0;ivalues[i]+=scalar*basis[i]; - } - } - delete gauss; - - /* Add residual if necessary*/ - gauss=new GaussTria(); - for(int iv=0;ivGaussVertex(iv); - - connectivity = IssmDouble(vertices[iv]->Connectivity()); - residual_input->GetInputValue(&residual,gauss); - pe->values[iv]+=residual/connectivity; - } - - /*Clean up and return*/ - xDelete(basis); - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorL2ProjectionEPL {{{*/ -ElementVector* Tria::CreatePVectorL2ProjectionEPL(void){ - - /*Intermediaries */ - int i,input_enum; - IssmDouble Jdet,value; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble slopes[2]; - Input* input = NULL; - Input* input2 = NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element vector*/ - ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&input_enum,InputToL2ProjectEnum); - switch(input_enum){ - case EplHeadSlopeXEnum: input2 = inputs->GetInput(EplHeadEnum); _assert_(input2); break; - case EplHeadSlopeYEnum: input2 = inputs->GetInput(EplHeadEnum); _assert_(input2); break; - default: input = inputs->GetInput(input_enum); - } - - /* Start looping on the number of gaussian points: */ - GaussTria* gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - - if(input2) input2->GetInputDerivativeValue(&slopes[0],&xyz_list[0][0],gauss); - switch(input_enum){ - case EplHeadSlopeXEnum: value = slopes[0]; break; - case EplHeadSlopeYEnum: value = slopes[1]; break; - default: input->GetInputValue(&value,gauss); - } - - for(i=0;ivalues[i]+=Jdet*gauss->weight*value*basis[i]; - } - - /*Clean up and return*/ - xDelete(basis); - delete gauss; - return pe; } /*}}}*/ @@ -7616,238 +6490,4 @@ } /*}}}*/ -/*FUNCTION Tria::CreatePVectorMasstransport{{{*/ -ElementVector* Tria::CreatePVectorMasstransport(void){ - - switch(GetElementType()){ - case P1Enum: case P2Enum: - return CreatePVectorMasstransport_CG(); - case P1DGEnum: - return CreatePVectorMasstransport_DG(); - default: - _error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet"); - } -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorMasstransport_CG {{{*/ -ElementVector* Tria::CreatePVectorMasstransport_CG(void){ - - /*Intermediaries */ - IssmDouble Jdet,dt; - IssmDouble ms,mb,mb_correction,thickness; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element vector and other vectors*/ - ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - Input* ms_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(ms_input); - Input* mb_input = inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(mb_input); - Input* mb_correction_input = inputs->GetInput(BasalforcingsMeltingRateCorrectionEnum); - Input* thickness_input = inputs->GetInput(ThicknessEnum); _assert_(thickness_input); - - /*Initialize mb_correction to 0, do not forget!:*/ - /* Start looping on the number of gaussian points: */ - GaussTria* gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - - ms_input->GetInputValue(&ms,gauss); - mb_input->GetInputValue(&mb,gauss); - thickness_input->GetInputValue(&thickness,gauss); - if(mb_correction_input) - mb_correction_input->GetInputValue(&mb_correction,gauss); - else - mb_correction=0.; - - for(int i=0;ivalues[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb-mb_correction))*basis[i]; - } - - /*Clean up and return*/ - xDelete(basis); - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorMasstransport_DG {{{*/ -ElementVector* Tria::CreatePVectorMasstransport_DG(void){ - - /*Intermediaries */ - IssmDouble Jdet,dt; - IssmDouble ms,mb,thickness; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element vector and other vectors*/ - ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* ms_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(ms_input); - Input* mb_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(mb_input); - Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); - - /* Start looping on the number of gaussian points: */ - GaussTria* gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - - ms_input->GetInputValue(&ms,gauss); - mb_input->GetInputValue(&mb,gauss); - thickness_input->GetInputValue(&thickness,gauss); - - for(int i=0;ivalues[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*basis[i]; - } - - /*Clean up and return*/ - xDelete(basis); - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorFreeSurfaceTop {{{*/ -ElementVector* Tria::CreatePVectorFreeSurfaceTop(void){ - - /*Intermediaries */ - IssmDouble Jdet,dt; - IssmDouble ms,surface,vz; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element vector and other vectors*/ - ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - Input* vz_input = inputs->GetInput(VzEnum); _assert_(vz_input); - Input* ms_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(ms_input); - Input* surface_input= inputs->GetInput(SurfaceEnum); _assert_(surface_input); - - /*Initialize mb_correction to 0, do not forget!:*/ - /* Start looping on the number of gaussian points: */ - GaussTria* gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - - vz_input->GetInputValue(&vz,gauss); - ms_input->GetInputValue(&ms,gauss); - surface_input->GetInputValue(&surface,gauss); - - for(int i=0;ivalues[i]+=Jdet*gauss->weight*(surface + dt*ms + dt*vz)*basis[i]; - } - - /*Clean up and return*/ - xDelete(basis); - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorFreeSurfaceTop1D {{{*/ -ElementVector* Tria::CreatePVectorFreeSurfaceTop1D(void){ - - if(!HasEdgeOnSurface()) return NULL; - - int index1,index2; - this->EdgeOnSurfaceIndices(&index1,&index2); - - Seg* seg=(Seg*)SpawnSeg(index1,index2); - ElementVector* pe=seg->CreatePVectorFreeSurfaceTop(); - delete seg->material; delete seg; - - /*clean up and return*/ - return pe; -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorFreeSurfaceBase {{{*/ -ElementVector* Tria::CreatePVectorFreeSurfaceBase(void){ - - /*Intermediaries */ - IssmDouble Jdet,dt; - IssmDouble mb,mb_correction,bed,vz; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element vector and other vectors*/ - ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - Input* vz_input = inputs->GetInput(VzEnum); _assert_(vz_input); - Input* mb_input = inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(mb_input); - Input* mb_correction_input = inputs->GetInput(BasalforcingsMeltingRateCorrectionEnum); - Input* bed_input = inputs->GetInput(BedEnum); _assert_(bed_input); - - /*Initialize mb_correction to 0, do not forget!:*/ - /* Start looping on the number of gaussian points: */ - GaussTria* gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - - vz_input->GetInputValue(&vz,gauss); - mb_input->GetInputValue(&mb,gauss); - bed_input->GetInputValue(&bed,gauss); - if(mb_correction_input) - mb_correction_input->GetInputValue(&mb_correction,gauss); - else - mb_correction=0.; - - for(int i=0;ivalues[i]+=Jdet*gauss->weight*(bed+dt*(mb-mb_correction) + dt*vz)*basis[i]; - } - - /*Clean up and return*/ - xDelete(basis); - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorFreeSurfaceBase1D {{{*/ -ElementVector* Tria::CreatePVectorFreeSurfaceBase1D(void){ - - if(!HasEdgeOnBed()) return NULL; - - int index1,index2; - this->EdgeOnBedIndices(&index1,&index2); - - Seg* seg=(Seg*)SpawnSeg(index1,index2); - ElementVector* pe=seg->CreatePVectorFreeSurfaceBase(); - delete seg->material; delete seg; - - /*clean up and return*/ - return pe; -} -/*}}}*/ #endif @@ -7982,64 +6622,4 @@ delete gauss; return Ke; -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorDamageEvolution{{{*/ -ElementVector* Tria::CreatePVectorDamageEvolution(void){ - - switch(GetElementType()){ - case P1Enum: case P2Enum: - return CreatePVectorDamageEvolution_CG(); - case P1DGEnum: - _error_("DG not implemented yet"); - default: - _error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet"); - } -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorDamageEvolution_CG {{{*/ -ElementVector* Tria::CreatePVectorDamageEvolution_CG(void){ - - /*Intermediaries */ - IssmDouble Jdet ,dt; - IssmDouble f,damage; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble f_list[NUMVERTICES]; - Input *damage_input = NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element vector and other vectors*/ - ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - damage_input = this->material->inputs->GetInput(DamageDbarEnum); _assert_(damage_input); - - /*retrieve damage evolution forcing function: */ - this->DamageEvolutionF(&f_list[0]); - - /*Initialize forcing function f to 0, do not forget!:*/ - /* Start looping on the number of gaussian points: */ - GaussTria* gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - - TriaRef::GetInputValue(&f,f_list,gauss); - damage_input->GetInputValue(&damage,gauss); - - for(int i=0;ivalues[i]+=Jdet*gauss->weight*(damage+dt*f)*basis[i]; - } - - /*Clean up and return*/ - xDelete(basis); - delete gauss; - return pe; } /*}}}*/ @@ -8550,284 +7130,4 @@ xDelete(basis); return Ke; -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorBalancethickness{{{*/ -ElementVector* Tria::CreatePVectorBalancethickness(void){ - - switch(GetElementType()){ - case P1Enum: - return CreatePVectorBalancethickness_CG(); - break; - case P1DGEnum: - return CreatePVectorBalancethickness_DG(); - default: - _error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet"); - } -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorBalancethickness_CG{{{*/ -ElementVector* Tria::CreatePVectorBalancethickness_CG(void){ - - /*Intermediaries */ - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble dhdt_g,mb_g,ms_g,Jdet; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element vector and other vectors*/ - ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* ms_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(ms_input); - Input* mb_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(mb_input); - Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input); - - /* Start looping on the number of gaussian points: */ - GaussTria* gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - ms_input->GetInputValue(&ms_g,gauss); - mb_input->GetInputValue(&mb_g,gauss); - dhdt_input->GetInputValue(&dhdt_g,gauss); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - - for(int i=0;ivalues[i]+=Jdet*gauss->weight*(ms_g-mb_g-dhdt_g)*basis[i]; - } - - /*Clean up and return*/ - xDelete(basis); - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorBalancethickness_DG {{{*/ -ElementVector* Tria::CreatePVectorBalancethickness_DG(void){ - - /*Intermediaries */ - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble mb_g,ms_g,dhdt_g,Jdet; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element vector and other vectors*/ - ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* ms_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(ms_input); - Input* mb_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(mb_input); - Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum);_assert_(dhdt_input); - - /* Start looping on the number of gaussian points: */ - GaussTria* gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - ms_input->GetInputValue(&ms_g,gauss); - mb_input->GetInputValue(&mb_g,gauss); - dhdt_input->GetInputValue(&dhdt_g,gauss); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - - for(int i=0;ivalues[i]+=Jdet*gauss->weight*(ms_g-mb_g-dhdt_g)*basis[i]; - } - - /*Clean up and return*/ - xDelete(basis); - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorBalancevelocity{{{*/ -ElementVector* Tria::CreatePVectorBalancevelocity(void){ - - /*Intermediaries */ - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble dhdt_g,mb_g,ms_g,Jdet; - IssmDouble h,gamma,thickness; - IssmDouble hnx,hny,dhnx[2],dhny[2]; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element vector and other vectors*/ - ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); - IssmDouble* basis = xNew(numnodes); - IssmDouble* dbasis = xNew(numnodes*2); - IssmDouble* H = xNew(numnodes); - IssmDouble* Nx = xNew(numnodes); - IssmDouble* Ny = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* ms_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(ms_input); - Input* mb_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(mb_input); - Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input); - Input* H_input=inputs->GetInput(ThicknessEnum); _assert_(H_input); - h=sqrt(2.*this->GetArea()); - - /*Get vector N for all nodes*/ - GetInputListOnNodes(Nx,SurfaceSlopeXEnum); - GetInputListOnNodes(Ny,SurfaceSlopeYEnum); - GetInputListOnNodes(H,ThicknessEnum); - for(int i=0;ibegin();igend();ig++){ - - gauss->GaussPoint(ig); - - ms_input->GetInputValue(&ms_g,gauss); - mb_input->GetInputValue(&mb_g,gauss); - dhdt_input->GetInputValue(&dhdt_g,gauss); - H_input->GetInputValue(&thickness,gauss); - if(thickness<50.) thickness=50.; - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss); - - TriaRef::GetInputDerivativeValue(&dhnx[0],Nx,&xyz_list[0][0],gauss); - TriaRef::GetInputDerivativeValue(&dhny[0],Ny,&xyz_list[0][0],gauss); - TriaRef::GetInputValue(&hnx,Nx,gauss); - TriaRef::GetInputValue(&hny,Ny,gauss); - - gamma=h/(2.*thickness+1.e-10); - - for(int i=0;ivalues[i]+=Jdet*gauss->weight*(ms_g-mb_g-dhdt_g)*( - basis[i] + gamma*(basis[i]*(dhnx[0]+dhny[1])+hnx*dbasis[0*numnodes+i] + hny*dbasis[1*numnodes+i]) - ); - } - } - - /*Clean up and return*/ - xDelete(basis); - xDelete(dbasis); - xDelete(H); - xDelete(Nx); - xDelete(Ny); - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorSmoothedSlopeX{{{*/ -ElementVector* Tria::CreatePVectorSmoothedSlopeX(void){ - - /*Intermediaries */ - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble Jdet; - IssmDouble thickness,slope[2]; - IssmDouble taud_x,norms,normv,vx,vy; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element vector and other vectors*/ - ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* H_input = inputs->GetInput(ThicknessEnum); _assert_(H_input); - Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input); - Input* vx_input = inputs->GetInput(VxEnum); - Input* vy_input = inputs->GetInput(VyEnum); - - /* Start looping on the number of gaussian points: */ - GaussTria* gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - H_input->GetInputValue(&thickness,gauss); - surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss); - if(vx_input && vy_input){ - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - norms = sqrt(slope[0]*slope[0]+slope[1]*slope[1]+1.e-10); - normv = sqrt(vx*vx + vy*vy); - if(normv>15./(365.*24.*3600.)) slope[0] = -vx/normv*norms; - } - taud_x = matpar->GetRhoIce()*matpar->GetG()*thickness*slope[0]; - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - - for(int i=0;ivalues[i]+=Jdet*gauss->weight*taud_x*basis[i]; - } - - /*Clean up and return*/ - xDelete(basis); - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorSmoothedSlopeY{{{*/ -ElementVector* Tria::CreatePVectorSmoothedSlopeY(void){ - - /*Intermediaries */ - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble Jdet; - IssmDouble thickness,slope[2]; - IssmDouble taud_y,norms,normv,vx,vy; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element vector and other vectors*/ - ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* H_input = inputs->GetInput(ThicknessEnum); _assert_(H_input); - Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input); - Input* vx_input = inputs->GetInput(VxEnum); - Input* vy_input = inputs->GetInput(VyEnum); - - /* Start looping on the number of gaussian points: */ - GaussTria* gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - H_input->GetInputValue(&thickness,gauss); - surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss); - if(vx_input && vy_input){ - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - norms = sqrt(slope[0]*slope[0]+slope[1]*slope[1]+1.e-10); - normv = sqrt(vx*vx + vy*vy); - if(normv>15./(365.*24.*3600.)) slope[1] = -vy/normv*norms; - } - taud_y = matpar->GetRhoIce()*matpar->GetG()*thickness*slope[1]; - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - - for(int i=0;ivalues[i]+=Jdet*gauss->weight*taud_y*basis[i]; - } - - /*Clean up and return*/ - xDelete(basis); - delete gauss; - return pe; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16981) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16982) @@ -222,20 +222,4 @@ ElementMatrix* CreateMassMatrix(void); ElementMatrix* CreateBasalMassMatrix(void); - ElementVector* CreatePVector(void); - ElementVector* CreatePVectorBalancethickness(void); - ElementVector* CreatePVectorBalancethickness_DG(void); - ElementVector* CreatePVectorBalancethickness_CG(void); - ElementVector* CreatePVectorBalancevelocity(void); - ElementVector* CreatePVectorSmoothedSlopeX(void); - ElementVector* CreatePVectorSmoothedSlopeY(void); - ElementVector* CreatePVectorMasstransport(void); - ElementVector* CreatePVectorMasstransport_CG(void); - ElementVector* CreatePVectorMasstransport_DG(void); - ElementVector* CreatePVectorFreeSurfaceTop(void); - ElementVector* CreatePVectorFreeSurfaceTop1D(void); - ElementVector* CreatePVectorFreeSurfaceBase(void); - ElementVector* CreatePVectorFreeSurfaceBase1D(void); - ElementVector* CreatePVectorL2Projection(void); - ElementVector* CreatePVectorL2ProjectionBase(void); IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");}; IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");}; @@ -306,13 +290,5 @@ ElementMatrix* CreateKMatrixStressbalanceFSViscous(void); ElementMatrix* CreateKMatrixStressbalanceFSFriction(void); - ElementVector* CreatePVectorStressbalanceSSA(void); - ElementVector* CreatePVectorStressbalanceSSADrivingStress(void); - ElementVector* CreatePVectorStressbalanceSSAFront(void); - ElementVector* CreatePVectorStressbalanceSIA(void); - ElementVector* CreatePVectorStressbalanceFS(void); - ElementVector* CreatePVectorStressbalanceFSFront(void); - ElementVector* CreatePVectorStressbalanceFSViscous(void); void PVectorGLSstabilization(ElementVector* pe); - ElementVector* CreatePVectorStressbalanceFSShelf(void); ElementMatrix* CreateJacobianStressbalanceSSA(void); IssmDouble GetYcoord(GaussTria* gauss); @@ -322,6 +298,4 @@ ElementMatrix* CreateKMatrixAdjointBalancethickness(void); ElementMatrix* CreateKMatrixAdjointSSA(void); - ElementVector* CreatePVectorAdjointHoriz(void); - ElementVector* CreatePVectorAdjointBalancethickness(void); #endif @@ -339,8 +313,4 @@ ElementMatrix* CreateKMatrixHydrologyDCEfficient(void); ElementMatrix* CreateEPLDomainMassMatrix(void); - ElementVector* CreatePVectorHydrologyShreve(void); - ElementVector* CreatePVectorHydrologyDCInefficient(void); - ElementVector* CreatePVectorHydrologyDCEfficient(void); - ElementVector* CreatePVectorL2ProjectionEPL(void); void CreateHydrologyWaterVelocityInput(void); void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index); @@ -356,6 +326,4 @@ ElementMatrix* CreateKMatrixDamageEvolution(void); ElementMatrix* CreateKMatrixDamageEvolution_CG(void); - ElementVector* CreatePVectorDamageEvolution(void); - ElementVector* CreatePVectorDamageEvolution_CG(void); void DamageEvolutionF(IssmDouble* flist); #endif Index: /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp (revision 16981) +++ /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp (revision 16982) @@ -42,5 +42,6 @@ element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); ElementMatrix* Ke = element->CreateKMatrix(); - ElementVector* pe = element->CreatePVector(); + //ElementVector* pe = element->CreatePVector(); + ElementVector* pe = analysis->CreatePVector(element); element->ReduceMatrices(Ke,pe); if(Ke) Ke->AddToGlobal(Kff_temp,NULL); @@ -72,8 +73,8 @@ for (i=0;ielements->Size();i++){ element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); - //ElementVector* pe = analysis->CreatePVector(element); + ElementMatrix* Ke = element->CreateKMatrix(); //ElementMatrix* Ke = analysis->CreateKMatrix(element); - ElementVector* pe = element->CreatePVector(); - ElementMatrix* Ke = element->CreateKMatrix(); + //ElementVector* pe = element->CreatePVector(); + ElementVector* pe = analysis->CreatePVector(element); element->ReduceMatrices(Ke,pe); if(Ke) Ke->AddToGlobal(Kff,Kfs);